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NOTE ON THE MOTION OF AN INFINITE CYLINDER 
. IN ROTATING VISCOUS LIQUID 
By W. R. DEAN (University College, London) 
[Received 27 November 1952] 


SUMMARY 
It has been shown by Sir Geoffrey Taylor that in certain conditions the two- 
dimensional motion of an infinite rigid cylinder in inviscid liquid is not altered if a 
uniform rotation is imposed on the whole system; it is shown here that the same 
result holds if the liquid is viscous. 


Sm GEOFFREY TAYLOR (1) has considered the two-dimensional motion of 
an infinite rigid cylinder in rotating inviscid liquid which is assumed homo- 
geneous and incompressible. Two motions were compared, the second being 
derived from the first by the superposition on the whole system of a rotation 
with constant angular velocity about a fixed axis, and Taylor has reduced 
toa very simple form the difference between the reactions of the liquid on 
the cylinder in the two cases. One conclusion, which has been confirmed 
by some striking experiments (2), is that, if the density of the cylinder is 
equal to that of the liquid, the movement under given forces of the cylinder 
relative to the liquid is not altered by the superposition of a uniform 
rotation on the whole system. 

Taylor’s argument is given in sections 1-3; in section 4 it is shown to 
apply with hardly any change to the motion of a cylinder in viscous 
liquid. 


1. In the first motion the pressure is denoted by p,, and (w,, v,) and (f;,9;) 
are the components, referred to rectangular Cartesian coordinates (x, y), 
of the velocity and acceleration of the liquid. Then the equations of motion 
of liquid of constant density p are 

phi ~Op,/OX, = py = —Op,/ey, 
and the equation of continuity is satisfied if we write 
U, = —Op,/doy, Vv, = Of, /dx. 
The condition at any point P of the boundary is satisfied if 


(Vin = —ap,/és, (3) 


where s, n are the directions shown in the figure and (V,), is the component 
along n of the velocity of the point P of the cylinder. 


{Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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The second motion of the liquid and the cylinder is derived from the 
first by superposing on the whole system a rotation about O with constant 
angular velocity w. The same coordinates (x,y) can be used to describe 
the second motion, but the axis Ox now makes an angle wt with a fixed 
axis OX. 











The particle of the liquid that is at the point (x, y) at time ¢ in the first 
motion is at the same point (a, y) referred to the rotating axes in the second 
motion. In the first case the velocity of the particle is (2, #/) and its accelera- 
tion (x, #/); in the second case the velocity is (¢—wy, y+ wx) and the accelera- 
tion (*¢—2wy—w*x, i+ 2wa—w*y). These are Lagrangian components; 
between the Eulerian components there are accordingly the relations 


Oy, 
Us = Uy —wy = —— — wy, (4) 
oy 
ous Pe 
Vo = Vi, +wr = - 1 se wer, (9) 
Cx 
: ‘ C P 
fe = fp —2wy—w?x = ite —w* nr, (6) 
Ox 
1 9,0 2 9 é 1 2 ‘| 
Jo = Jy 2wkL—wry = g,—2w— — wy. (/) 
ou 


The equation of continuity is satisfied (as it evidently must be) in the 
second motion, for 


Us = —Oy,/Cy, V> = Oy,/éx, 
where hy = y+ dw(x?+y?). (8) 
From (1) and the corresponding equations with suffix 2 we have 


O(P2—p)/Cx = —p(fe—f), APs—P,)/Cy = —p(J2—91): 
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MOTION OF AN INFINITE CYLINDER 


and then, using (6) and (7), 


( af; 
ri lo ),) Ow ° c 9/9 9 
Po P13 p| 2w 14 wr} = p— | Zw, + Ley? (x +y*) |, 
Cx CX Cx ” 
(Ps P) . ous, 9 c ) 1.2/2} 22 
= p| 2w + wy) = p— | 2wib, + bw?(a?+-y?) |; 
cy cy cy is 


these equations give Taylor’s formula for the pressure in the second motion, 
[rs 2/421 92 
Po = Py+ pl 2a, + bw?(x?+-y")]. (9) 


2. Taylor next compares the reactions on the cylinder in the two motions. 
Suppose that (/), G,) are the components of the force exerted on unit length 
of the cylinder in the first motion and that H, is the moment of the liquid 
pressures about C, the mass-centre of the cross-section. In the second 
motion (43, G,, H,) differ from (F,, G,, H,) because of the difference between 
the pressures p,, >. Hence, if L denotes the boundary of the section of the 


cylinder 


F,—F, ( (Po—p,) dy = —2pw i ys, dy — S pw" (a?+-y?) dy, 
L (L) (L) 


Gy—G, = | (Po—p,) dx = 2pw | f, dx + dpw* | (x+y?) dz. 
i L) (L) 


To evaluate the first line-integral Taylor writes 


j , r , oy ° Gus 
ys, dy o,— ds [by y|— | y— ds. (10) 
. + os . cs 
I I L) 
as . 
. CW 7 
From equation (3), | —*ds = - (V3), ds, 
Cs 
l (L) 


and the integral on the right is zero since the cylinder is rigid; thus there 
is no change in ys, in going once round L, and the integrated term in (10) 
shown by the bracket is zero. 

The cylinder is in two-dimensional motion; at time ¢ in the first motion 
let 2 be its angular velocity and let (a, yy) be the coordinates of C. Then 
at iny point P of the boundary 

Cus 


l V3) (%y9—Qy)cos x + (Yo + QE)sin x: (11) 


Cs 


here y is the angle between » and Ox, and (€, 7) are the coordinates of P 


relative to C, so that 2 y+, y Yo 9- 


U 
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From (10) and (11), 


cs 


[ vdy=— [ yHras 


(L) (L) 


| (Yo+7)(%—Qy) dy — (Yo n)(Jo+QE) dé 


7 | E (Yo + 9)(%y—Qy)} +5 uot mio + OE }| dédn, 


A is the area of the cross-section. Hence 
[ by dy = = | { ¢ Yor Q€) dédyn = = Ajo, 
(L) e 1) 
and similarly Taylor finds 
| #, dx = Axy. 
(L) 


The evaluation of the other line-integrals is straightforward: 
& =) 


| (x?+-y?) dy = | | 7 (x?+-y?) drdy = 2A2,, 
Ca 
(L) (A) 
| (x?+y?) dx = — | 2 (x?+-y?) dxdy = —2Ayp, 
A cy 
(L) (A) 
so that finally F.—F, = —pA(eo%ty-+ 2), (12) 
G,—G, = —pA(w*y,—2wi,). (13) 
Taylor’s next result is H,—H, = 9. (14) 


For the difference between the moments we have 


H,—H, = ( (P2—P,)(n cos xy — Esin x) ds 


(L) 
=p | (dé + 4 dn)[2ay+$o%(a?+y")]. 
(L) 


It is easy to show that the term in w? is zero, for 


2 


. 


(L) 


(x?+-y*)(E dé + ny dn) = | | E n(x?+-y?)} — = fE (ax? +") dédn 
o€ ' on 
A) 


- fl [2n(2y)+€)—2€(yp+n)] dédy 
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MOTION OF AN INFINITE CYLINDER 
In evaluating the other term Taylor first establishes the equation 
. r “ - 
W,(2€ dé + 2n dyn) = — | (2442) ds, (15) 
cs 
L) (L) 
which corresponds to (10), the integrated term again being zero, and then 
uses the boundary condition in the form (11). This gives 


[ y,(2 dé + 29 dn) = | (€2+y%)[(@>—Qy) dq— (yo ME) dé] 
i L) 


uy 


+n) (do+08)} dédy 


and equation (14) follows. 


3. Equations (12), (13), and (14) show that the forces due to the pressure 
on the cylinder in the second motion are the resultant of: (i) the forces that 
uct on the cylinder in the first motion, and (ii) a force 

(—pw*Ax,—2pwAiy, —pw*Ayy+2pwAxy) 
acting at C. 

Lastly Taylor considers the related dynamical problem: a uniform rigid 
cylinder of mass VW per unit length is made to move in any way in two dimen- 
sions, and then made to move in the same way relative to rotating axes. 
The forces required in the second movement are statically equivalent to: 
i) the forces required in the first movement, and (ii) a force 

(—Mw*2,—2Mwy,, —Mw*y,+2Mwiy) 


acting at C, C being defined as before. 

Hence if a cylinder immersed in liquid of the same density is moved in 
two dimensions by given external forces, then a uniform rotation of the 
whole system (including the external forces) will not alter the motion of 
the cylinder relative to the liquid. 

4. Suppcse now that the liquid is viscous, but homogeneous and incom- 
pressible as before. In this case the equations for the first motion are 


,' P ; - 
; CP, o* o* CP, , - . = 
t , 
Pji + | — = Uy, (, =H — mw =.” V1; 
oa | aq? aa . pg cy a éy?} * 
where p, is now the mean pressure. Equations (2) and (3) need no change, 


though (3) does not now give the only condition that must be satisfied at 
the boundary. If, however, there is no velocity of slip in the first motion, 
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there will be none in the second. Equations (4) to (8) are kinematical and 
still apply. Lastly, since 


9 


j 92 e2 \ C 2 e2 
= i) v1) = te ses Vy), 
Ox OY” a ee 


O(Pps—Pp,)/ex = —p(f.—f,), O(Pp2—Pp,)/Cy = —p(Jo—J)), 


so that Taylor’s equation (9) is unchanged in form, though it now gives 


it is still true that 


the relation between the mean pressures. 

We have now to calculate, as in section 2, the reactions on the cylinder 
in the two motions. Let (p,.,),, (p,,)1, and (p,,,,), denote the stresses in the 
liquid in the first motion; then 


cu 


9 1 c Uy c Vv) ae. V; 
(Prr)i = — Pit = ’ (Pry) > B Ty ’ (Pyy)1 : Pit <2 ’ 
Cx ‘ cy Cx : cy 
(16) 
and the components of the force exerted on an element ds of L are 
U(Prr)1 COS X+(Pyy), SiN x, (P,,), COS X+(Py,), in x} ds (17) 


per unit length of the cylinder. Now it is easy to see from equations 
(4), (5), and (16) that 


(Prr)2—(Prx) (P2—P1) = (Pyy)o—(Pyy)v> (Pryde = (Pry)r: 
hence the difference between the two stress-systems arises solely from the 
difference between the mean pressures. Expression (17) and the corre- 
sponding expression with suffix 2 then show that there is no change in the 
expressions for (F,—F,, G,—G,, H,—H,) as line integrals. 

It remains to consider the evaluation of the three line-integrals which is 
given in full in section 2. There are the two transformations (10) and (15), in 
each of which the integrated term vanishes, since ys, does not change in a 
single circuit of the boundary; the boundary condition (11) is used; and the 
integrals round L are expressed as integrals over A by the formula 

Pdzx+ Qdy = | | (' @ om ; “| dxdy, 
‘ JJ \OU oy 
(L) 1) 
all the functions P, Q that are concerned being polynomials. It is clear 
that in none of these stages is any change needed if the viscosity of the 
liquid is allowed for, and it is concluded that Taylor’s result, quoted in 
section 3, holds for a viscous liquid. 


REFERENCES 
1. G. I. Taytor, Proc. Roy. Soc. A, 93 (1916), 99-113. 
a. ibid. A, 100 (1921), 114-21. 
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COMPRESSIBLE SUBSONIC FLOW IN TWO- 
DIMENSIONAL CHANNELS WITH MIXED 
BOUNDARY CONDITIONS 
By L. C. WOODS (R.N.Z.A.F. Scientific Defence Corps) 


Received 6 February 1953] 


SUMMARY 


Equations are given for the calculation of the subsonic compressible flow in two- 


limensional channelst when the boundary conditions are ‘simply-mixed’, i.e. when 


the pressures are known over some contiguoust sections of the walls or bounding 
treamlines, and the shape of the remaining sections is assigned. Solutions to simple 
problems of this type in incompressible flow are occasionally met in the literature 
n fluid dynamics (see (1) and (2)). The method of solution is almost invariably based 
nan application of the Schwarz—Christoffel theorem on the conformal mapping of 
, polygon, an approach which is limited to problems in which the boundary con- 
litions are specified as step functions. The general solution obtained in this paper 
includes as trivial examples many of the solutions obtained by separate and often 


engthy applications of the Schwarz—Christoffel theorem, and, more important, 


provides a method of dealing with compressible subsonic flow past continuously 
urved walls. In this latter case the solution is given by an integral equation, which 
an be solved by a rapidly mverging iterative process. Five examples are solved 
n section 5, namely 


i) the flow of a jet deflected round a given curved surface (the Coanda effect, 





see (3)): 
ii) the design of a bend in a channel; 
iii) the flow of a stream up a step, with boundary-layer separation ; 


iv) Helmholtz flow of a jet impinging normally on a flat plate ; 


a generalization of Borda’s mouthpiece. 


1. Introduction 
} Nome nelature 
r,y) the physical plane. 
v+iy, 0 \ 1). 
n,s distances measured normal to and along a streamline respectively. 
7,4) velocity vector in polar coordinates. 
P; py local and stagnation densities respectively. 
¢,4) plane of velocity equipotentials (f = constant), and streamlines 


constant), where 


dd qds, dis p q dn. (1) 


Po 
By hannels’ here is meant any simply connected region of flow bounded by two 
streamlines, one of which can be at infinity. 
+ If necessary, via the point at infinity. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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U velocity at an infinite distance upstream (f = —oo). 
M local Mach number. 
B = (1— M?)!. 
r is defined by the equation 
4 v 
r= | (] —M?}{log—) (2) 
’ q 
q=U 
m is defined by the equation 
m = (1—M?2)} Po, (3) 
p 
w is defined by the equation 
w= d+im,w. (4) 
ve is a suffix to denote values at infinity upstream. 
h channel ‘width’ in the w-plane. 
H channel width at infinity upstream in the (2, y)-plane, thus from 


(1), (3), and (4) h = HUB... (5) 

The independent variables of the theory are ¢ and ¢ as defined by (1). It 
is convenient to combine them into the complex number w defined by 
equation (4), which in the case of incompressible flow becomes the well- 
known complex potential function. Since the boundary conditions in a 
given problem are usually specified, not in the w-plane, but in the z-plane 
of physical coordinates x and y, the use of ¢ and ¢ as independent variables 
leads in general to integral equations for the dependent variables (see, for 
example, equation (35)). However, integral equations appear to be un- 
avoidable in any exact theory dealing with curved walls. 

The dependent variables of the theory will be r and 0, where r is defined 
by equation (2) and @ is the direction of flow. It is shown in the Appendix 
that when the approximation 


m= mm, (6) 
is admissible, r and @ are conjugate harmonic functions in the w-plane.t 
This theory is only applicable to subsonic flow (see equation (2)), but is 
reasonably accurate provided M is less than the critical value corresponding 


to the first appearance of sonic speed locally. In compressible flow, then, 
the complex number defined by 


+ Equations (2) and (6) are due to von Karman, who, in (4), used them to show that ¢ 
and are approximately harmonic functions in the (r, 6)-plane. 


is Op) 


r=] 


Sinee 
follo 
flow. 
W 
inco! 
dedi 
how 


whi 
beti 
in ¢ 


sen 


na 








from 


(0) 


= 
l by 
vell- 
in a 
lane 
bles 
for 


un- 


ned 
dix 
(6) 


le. t 


ing 


en, 


(7) 


it d 














COMPRESSIBLE SUBSONIC FLOW IN TWO-DIMENSIONAL CHANNELS 265 


is approximately an analytic function of w. If the flow is incompressible, 
r= log(U/q), and 


, U 6 dz 
oi - pe = ] - . 8 
J . | q ' ) og(t i: (8) 


Since the analytic character of f is the basis of the theory of this paper, it 
follows from equation (8) that the theory will be exact for incompressible 
flow. 

When f(w) has been calculated by the methods given below in the case of 
incompressible flow, equation (8) enables z = z(w), and hence f = f(z), to be 
deduced, thus completing the solution. When the flow is compressible, 
however, the more complicated relation 


id 
dz = —(ag+ite ay), (9) 
q p 


which follows from equations (1), must be used to calculate the relation 
between the z- and w-planes. Usually, however, we shall only be interested 
in calculating the solution on the bounding streamlines, where 


? 
“1 
= _ dd. (10) 
J q 
The infinite strip —coo <d < 0, 0 < m,¢ < h, in the w-plane repre- 


sents flow in a channel of quite a general type. The bounding streamlines, 
w = 0and w = th, are lines on which either the pressure or the wall coordi- 
nates are known. With the aid of Bernoulli’s equation and equation (2), r can 
be calculated where the pressures are given, while 6 (6 = tan-1(dy/dz)) 
follows directly in the sections where the wall coordinates are given. Fig. 1 
shows a number of possible w-planes—the full lines showing intervals in 
which @ is known, and the dotted lines intervals in which r is known. The 
case h = oo is only a special case of channel flow. Now Figs. 1(a), 1(5),..., 1(g) 
all have the property in common that the sections of the walls on which 6 
or r) is known are contiguous—the point at infinity being a point on each 
wall. The method of this paper can only be applied to w-planes of this type. 
Fig. 1(h) is an example of a w-planet outside the scope of the method of the 
present paper—w-planes of this type will be considered in another paper. 
Figs. 1(a) and 1(g) show w-planes in which the boundary conditions are not 
mixed. These latter cases (the ‘indirect’ and ‘direct’ problems respectively) 
have been treated at length in (5). 

In the next section an important auxiliary plane ¢ = 5+ 7€ is introduced. 
This plane is a conformal transformation of the w-plane such that the 


+ A simple problem with a w-plane of this type is solved in (1). 
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intervals of the bounding streamlines on which @ is known are mapped on 
to € = 0, —%2 <8 < @, while the remaining intervals are mapped on to 
€ = in/2, —0 <8 < ~w. The ¢-plane will thus have the appearance of the 


w-plane shown in Fig. 1(f). The mapping equation ¢ = ¢(w) can be easily 














ke es ee: ee eee 
es... Cee oes See 
(a) (b) 
eee ATT ae Se 
ene a ee 
(c) (d) 
(e) (F) 
(a) (hy 


calculated for the cases (5), (c),..., (f) in Fig. 1 by making use of the Schwarz- 
Christoffel theorem. Then, since f(w) is analytic, so also is f(Z), and hence 
in bo anes - 
in both planes Vf = 0. (11) 
The solution of equation (11) in the comparatively simple ¢-plane is obtained 
in section 3. When this has been found, f = f(w) follows immediately. In 
the case of incompressible flow it follows from (8) that 


Uz = ( expf(w) dw, (12) 


but the integral can be calculated analytically only for the simplest prob- 
lems. For curved boundaries the solution of equation (11) is an integral 
equation, and so usually numerical methods have to be introduced at this 
point, f(w) then being obtained as a discrete set of numerical values. 


2. The auxiliary ¢-plane 
Fig. 2 illustrates the relation between the w- and ¢-planes for the case 
shown in Fig. 1(b). By an application of the Schwarz—Christoffel theorem 


it is easily proved that the planes are related by the equation 


TW, Tw 
tanh ¢ = coth{| —2°)}tanh|— ). (13) 
. ( 2h ) (F] 


If wy in this equation is real and equal to ¢y, say, then we have the case 
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shown in Fig. 2, but if wy is made equal to ¢,+7h, then we will have a 
»-plane of the type implicit in Fig. 4. Itis unnecessary to consider separately 
the cases shown in Figs. 1 (c), (e), and (f), for when the solution f = f(w) 
has been found from f = /f(¢) by using equation (13), the solutions for these 
planes can be deduced as special cases as follows: 

Fig. l(c): Take wy = do, shift origin to w do, let by > 2. 

Fig. 1 (¢ Take w, d+ih, replace w by ¢,+ih—w, let dy > ©. 





Fig. 1(f): Take w, dy, let dy > OO. 
w piane 
| d6=0 
‘Be w=th 
Dex Cx 
ke by —> ¢; > 
“BR 0 A w-0 
Pepe & 7 
TGs by.) 
Bx rf 2h Rips 
C-0 
7] 
Fic. 2 


When the w-plane is of the type shown in Fig. 1 (d), the mapping equation 
sfound to be 
ad rh 
tanh ¢ sech| —_— tanh| "), (14) 
h Lh 


) 


24) is the potential difference between points A and B. 


where 


3. The solution of Laplace’s equation in the (-plane 


THEOREM. If y(Z) .+i8 is harmonic in the strip 
NY : 5 : 7). € €. 
there C 6+21€, then provide dv satisfie s conditions (a) and (b) give n below. 


yatany point P(C) within the strip is given by 


l | a " 7 , ~ ~ 
-osech (o* OC) 2 sech (o* c)| dd*, (15) 
2« | Ye Ze | 
vhere 8, lim 8(8*, €) and «, lim «(4*,e—&). The conditions referred to 
= —>() &—0 
a) lim } R, é)e—7hi2e 0, (16) 


39 and x. are continuous functions except possibly at a finite number of 
points 8... s a Aare m, where they may have either logarithmic singularities 


or simple di S¢ ontinuitie S. 
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At such a singularity, on € = 0 for example, 
y ~ K log(€—5,) = K log|€—5,|+7K arg(¢—8,), 


and so if K is real (imaginary), By («)) has a simple discontinuity of Kz at 
5 = 6,, while ap» (fy) has a logarithmic infinity at 6 = 6.,. 


Proof. Suppose y is defined outside the strip —00 <8 < wH,0< E<e 
by the equation (3, €-+2en) = (—1)"y(8, é), (17) 


where n is a positive or negative integer, then y is an analytic function 
within and on the rectangles 


—R<i< R, 2nme<E< (2n+1ye, (18) 


provided they are suitably indented to exclude the logarithmic singularities 
of y. It follows from Cauchy’s theorem when n + 0, and from Cauchy’s 
integral when n = 0, that at a point P(¢) within the rectangle 





2 
n=0,y(f) 1 [ (y(5*, 2ne+0) — y(d*, (2n+1) Je ds 
n#~0, 0 2m J |S*4+2nte—C 8*4(2n+1)ie—ZJ 
—R 
+ r { y(R, 2ne+é*) y(— R, 2ne+-€*) \ age 





} \R42Qnie+ig*—C —R+2nie+ié*—f 

0 
as it is a matter of simple analysis to show that as the radii of the indenta- 
tions tend to zero so do the contributions to the contour integrals from these 
indentations. By addition, and with the aid of equation (17), we find 





R xy 
LPS f(—1)"(8*, 40) (—1)"718*,e—0)| 5 
v0) = 5 ‘2 2. \ Spaniel FE An tI ye—z Ot 
—R X=- 





y \R+ Qnie+ié*—C —R+ Saete—1] 





. NV na * n 
=| 2 > jf (—1)% y(R, €*) (—1)"y(—R, €*) | ae. (19) 


From a result given on page 147 of ref. (6), it follows that 
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eae ,3 (= a) = ee ( b -*) 











COM 


and § 


Ifth 
an € 
witl 
side 
desc 


Ad 


we 


wl 
th 








K7 at 


uN 


nction 


(18) 


arities 


ichy’s 





4¢ 


lenta- | 
these 
id 


(20) 











COMPRESSIBLE SUBSONIC FLOW IN TWO-DIMENSIONAL CHANNELS 269 


and so, taking lim of equation (19), we find 


N-= 
R 
l j ' oe i 
y(C) 1 (5*, -(0)cosech a (o*—C)— 
ee die J | 2e 
R 
; i oe 
-y(8*, «—0)cosech — iit dd* +- 
= § " 7 . “ 
1 R, €*)cosech ~ (R+1€* —f)— 
te | Ze 


-y(— R, €*)cosech id (—R+2€*-— } dé*, (21) 
am€ 


Ifthe same argument is applied to the conjugate point P(Z), instead of P(Z), 
an equation similar to (21) is obtained, but with ¢ replaced by Z%, and 
with the left-hand side of the equation replaced by zero, since P(Z) lies out- 
side all the rectangles defined by (18). The conjugate to the equation just 


described is 


R 
| i? = T a» 
= - : y(o*, 1 ()cosech (o*—C)— 
fie J | Ze 
PR 
_ a ae ‘ 
y(5*, e—0)cosech (8*—ie—z)| dd* + 
Ze | 
l f fap ex vad 72 an 4 
y(R, €*)cosech — (R—1€* —C)— 
t¢ J | Ze 


y(— R, €*)cosech = (— R—ige—0)| dé*, (22) 


Adding equations (21) and (22), and making use of the results 


y(6*, +0)—y(d*, +0) 218, y(5*, «—0)+7(6*, «—0) = 2a,, 
7 . a a TT se y 
cosech — (6*+te—f) = Fisech — (5*—{), 
ué ae 


we obtain 


a (a, cosech (8*—{)+-a, sech  (8*—{)| d3* + I, 
2e 2e j 


where J is the contribution due to the integrals from 0 to «. To complete 
the proof of the theorem it only remains to let R tend to infinity, when, 
from condition (a) of the theorem, it follows that J tends to zero. 
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The conjugate equation to (15) is (replace 78) by ag and «, by if.) 


) 


pe 


x, cosech = (5*—Z)—B, sech hd (8*— 
| 2e 2e 


dd*. — (23) 


< 
> 
8 _——+, g 


We note in — that the equation 


~ ] ° 7 ~ 7 a -” 
ve) = 5 | {Po coth 5 (5*—0)—B, tanh F (5*—0) d3* +. 


+4{a(00)+a(—0o)}, (24) 
and its conjugate 


0 coth 5 (8- a xetanh F (5*—<()| d3* +. 


y(t) = —> 





bo 
ra) 


; 
' f ” * 
+5%B(%)+-B(—2)}, 
which have been used to solve the direct and indirect problems in ref. (5), 
‘an be proved in just the same way as the above theorem. The only differ- 
ence is that equation (17) is replaced by y(85,é+2en) = y(8, €), and instead 
of equation (20) it is necessary to use 


Ny 

: l 7a 
li + ( oth | — 
ed a. a cab) (Feo ( b } 


Equation (15) is essentially the solution of Laplace’s equation in an 
infinite strip when boundary values are given on € = 0, and normal deriva- 
tives are given on € = «. This can be shown as follows. Integration of 
equation (15) by parts gives 


x 
. 


9 - 
y(f) = a4 t+tB,+- tanh exp (*—1)} dp,(d*)— 


78 


tan exp = (3* 0) dx,(6*), (25) 
€ 





where a, and 8 , are the values of «and £8 at 6 = oo. When «, is continuous, 
the Cauchy—Riemann equations relating « and f yield 


da, = (<s) ase = (“F*) dd*. 
00 * o€ 


The statement made about equation (15) follows when this result is substi- 
tuted in equation (25). 
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4, The general solution 
The functions r and 6, defined in section 1, are conjugate harmonic 
functions corresponding to a and 8 of the theorem of section 3; thus, 
assuming conditions (a) and (b) of the theorem to hold, it follows from 

equations (7) and (15) that 
F(¢) . {6, cosech(5*— Z)+-r, sech(5* —)} dd*, (26) 


7i 
« 
x 


where, without loss of generality, « has been given the value 7/2 as in 
section 2. (However, the subscript ¢ will be retained for convenience of 
notation.) An alternative form of equation (26) is (cf. equation (25)) 


9 + 
f(Z) +10, + = tanh-lfexp(d*—)} d,(8*)— 
7 
2 7 m is 
_ tan-exp(d*—)} dr.(6*). (27) 
"he a 


When the w-plane is of one of the types (b), (c), or (f) shown in Fig. 1, 
then exp(6*—) may be calculated as follows. From equation (13), 








, sinh $a(wo+w) )/h\4 
expc 

\sinh x(w,—w)/hJ 

| 1 

sinh $n(w,+w*)/h 
and exp 6* —- > [2 

| sinh 4ar(wy—w*)/h & 
where w = w* when ¢ = 5*. Thus the term exp(5*— C) in equation (27) can 


be replaced by 
sinh $7(wy+w*)/I : |4 (sinh }27(wy—w)/h)} 


\sinh 4(wy+-w)/h} 


where w* is either d* or 6*+ th thers on the range of integration. 


exp(o*—C) 


(28) 











| sinh Lar(wy —w -w*) 


Similarly, as mentioned in section 2, w, is either ¢, or ¢,+th. The various 
special cases are now obtained as described in section 2. 
Similarly for the case shown in Fig. 1 (d), equation (14) yields 
exp(8*—Z) exp(—7¢o/h)+exp(mw*/h)|$( exp(7o/h)—exp(zw/h) \} 
; exp(7¢)/h)—exp(zw*/h) | \exp(—7do/h)+exp(aw/h)| 
(29) 


Discontinuities 7; in 6) and o; in r, at 8* = 8, plainly contribute terms 








> » 
y «0; y 
tanh-Yexp(5;—¢)} and ——/tan—{exp(5;—Z)} 


77 


to the right-hand side of equation (27). These terms can be expressed in 
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terms of ¢; (corresponding to ;) and w by means of equations (28) or (29). Dis- 
continuities in 6 (sharp corners are of course common in channel problems) 
lead to infinite values in the conjugate function r, but discontinuities in ; 
are really only of theoretical interest, since they lead to physically im- 
possible infinite spirals on the surface (see, for example, equation (49) in 
section 5). 

Where wall pressures, and hence, with the aid of Bernoulli’s theorem. 
velocities, are given as continuous functions of s, 





$(s) = U | 1 ds, (30) 


and equation (2) enables dr.(d*) to be determined and substituted in 
equation (27). However, where wall shapes are given, we can only write 
db,(4*) = So So qge  _ a" (31 
ds, d* Rod 
where R, is the radius of curvature, —(ds)/d6,). Thus in the region where 
the shape is given (—oo < 8 < o&) equation (27) yields 








po 
9 ¢ 
ro(¢) = ‘= tanh-l{exp(5*—Z)} ie 
ot=—do 0710 
9 J 
“4 i tan-fexp(o*—O)(775) dg*+- 


$*=—do 
+-terms due to discontinuities in #, and r,, (32) 


where and f denote integrals taken over the region in which are given 


the wall shape and pressure respectively. 

Equation (32) is an integral equation, q) appearing implicitly on the left- 
hand side and explicitly under the first integral on the right-hand side, and 
it must be solved in general by an iterative process. This process, which is 
illustrated by the first example in section 5, is as follows. Approximate 
values of g)(¢) are assumed and substituted in 

1 * U dd 

806) = FF | i (33) 
to yield an approximate relation s)(¢). Then, from this approximate 
relation and the given values of 4,(s), R(¢,), and hence R, qo, can be obtained 
approximately. This is substituted in equation (32) permitting r,(¢), and 
hence, from (2), g9(4), to be calculated. This completes the first iteration. 
This calculation is repeated until q,(¢) is unchanged from one iteration to 
the next (assuming convergence). 
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5. Examples 

In this section a variety of examples, each of some intrinsic interest, will 
be calculated to illustrate the application of the theory of the previous 
sections. 
5.1. The flow of a jet deflected round a given curved surface 

This example has a w-plane of type (b), Fig. 1. Writing wy = ¢o, we find 
from equations (27) and (28) that 


f(w) r,+20,-4 








A A! 
2 — ,{sinh{3a(¢9+¢*)/h} jsinh{}7( dy—w) VOH\3 a9 ($*) 
‘- J \sinh{3z(d, b*)/h isinh{42r(d9+ w)/h\| . 
_3| } } ltan ,/sinh{ a77(p* + ~$o)) hi jsinh{$77( go—w)/h}\4 * dro(¢*)— 
m st \sinh{42(6*—dy) )/hisinh{4a(d>+w)/h}} ’ 
o* g* 0 
2 - ees b77($o+$*)/hjsinh{}7(¢)—w) [AS\3 4, (*), (34) 
7 . \cosh{4z(d5 b*)/h} jsinh{$2(¢d9+w) hi} . 


where we have separated r.(4*) into ro(¢*), the value of r on w = 0, and 
r,(d*), the value of r on w = th. 


In the special case of a constant pressure jet, 7) = 7; = 0. In this case, 





with the aid of (31), equation (34) becomes on w = ¢, —d» < d < do, 
do 
(6) 2-7 |] — _(sinh{32(d9+¢4*)/h} isinh{3n($o bo- zA ta ” 
a Rod \sinh{32( (do d*) h' jsinh{4 \n(dot hi} 35) 


in which it has been assumed for simplicity that there are no sharp corners 
on the surface. In general, equation (35) must be solved numerically. 
A suitable scheme will now be described. 

The range (4), —¢,) is subdivided into (¢,,.4,6,—45-+-5 1454) 80 that 
(1/R,q)) can be taken as constant in each interval with negligible error; 


then we can write 


db, . 2~— ‘ 
W--3 2) eg 
ne ,(sinh{37(¢9+4*)/h} sinh{}7(o —;) M3 age 
\sinh{4(¢,—¢*)/h}sinh{42($9+-¢,) h\} , 
ie. r(¢;) = - = S (3) Au (36) 
| ¢ 








L. C. WOODS 


where (#/Rq); is the value of (49/ Ry qo) in the middle of the range (4, , +> 9;-1), 








$i+3/bo 
; ; tanh-! sna 2h}sinh{mdo(1— i 2h})4 it 
' ite \sinh{7¢,(1— )/2h} 'sinh{ado(1-+ 2h} 
i—3/Po 


and t; = ¢;/¢ 9. Suppose that the tangents at the points of first contact 
and separation of the jet include an angle 7 (Fig. 3), then 


$o 
d6(p*) = 7, 
ey 


1 ? 
and hence J ( ddé* = —r, 
$o Rq 


do 


i.e. 7 (Rt) ($i14—$i-4) +7 = 0. (37) 


5.11. ‘Exact’ numerical solution of the incompressible flow of a jet past a 
circular cylinder 


As a particular example, consider the flow of an incompressible fluid past 


a circular cylinder with + = —7/3, and the jet width such that 
_ Tho 29 
= —f=— 1, (38 
2h 


An appropriate matrix A,;, which was calculated for the author by the 
Mathematics Division of the National Physical Laboratory, is given in 
Table 1. The range (¢),—¢ 9) has been divided into ten equal intervals. 
From equations (2), (5), (36), (37), and (38) we find 
ul 
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q : 10 ¥ (U/q), A; 
Thus logs (d;) > are __. (40) 
d 2 
- m2 (Y/N: 

ABLE | 

By A.. 

108 x A;; 

5 6 7 } 5 9 10 
2 I 7 357 296 224 166 113 58 
1178 17 172 r2¢ 819 614 462 339 231 11g 
3502 4¢ 1288 933 689 501 340 174 
1129 1947 54 2045 1353 965 689 462 | 236 
sot 82 1288 048 3716 2078 1353 933 616 | 310 
10 61 933 5 2078 3716 2048 1288 820 | 406 
23¢ 462 68 65 1353 2048 3654 1947 1129 539 
17 34 ol 933 1288 1946 3502 1723 747 
I 462 614 819 1126 1722 3179 1175 
113 11 4 296 387 516 716 1138 | 2200 
‘ ‘ re 





Toobtain a first approximation it is assumed that g = U, when (40) becomes 


log 4 lS 4. 
eit 2 Bo 


and with the aid of Table | we arrive at the first column of Table 2. (The 
values of g/U at i = 1...., 5 follow by symmetry.) If the values of q/U in 
this column are used in the right-hand side of equation (40), there result 
the values given in the second column of Table 2, and so on. The con- 
vergence is clearly quite rapid. Column 5 of Table 2 gives the angular 


TABLE 2 


Velocity distribution 





2 3 4 | 5 6 

q/l q/U w q/l 
5 I 1*534 1°539 1*535 | 2°76 1*570° 
¢ 1°542 "513 1°516 516 | 8&2 1°547 
7 1°482 4 1*468 1°468 | 14°04 1°494 
8 1°384 84 1°384 1°384 20°00 1*404 
I 31 1*230 1°230 26°54 1°237 
Ic 1*00C 000 1*°000 1°000 30°00 1°000 





coordinate measured from the midpoint of the ‘wetted’ circular arc. It was 
obtained by use of equation (10), Viz. 
d 


sU 1fU ~ LX [U) (bin: 
2¢ 5 | q ae) 73 (2). (2 x =) - 5: 


0 / 


From column 4 of Table 2 and equation (39) it is found that the depth of 


the jet at infinity is given by 


H 1-166R. 
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5.12. Approximate analytical solution of the flow of a jet past a circular 
cylinder 
If g under the integral sign in equation (35) is replaced by ¢, the mean 
value of q in (¢5, —4 9), then in the ¢-plane the equation becomes 
x 


a 


r(d) = : , tanh-lfexp(s*—8)! a" 15* 
. ear € 2x - ao*. 
, aRq eX] Vase] 


© 


With the aid of equation (13) we find 
> 2hl 
. re) a eee as 
re) a Rg 
where ‘a 
_—_ | tanh-{exp(5*—8)} d8* 


cosh25*— tanh?(7¢,/2h)sinh?* 











7 coth To tan Yo a 

2h | coshd = }’ 

‘ 2h (sinh(29/2h)) 

reyes 48) <2 tant Pol =M)\ 
whence r(d) = Ri un oe oe 


If s) is the surface distance on the cylinder from ¢ = 0 to d = dp, then 
98 = 7'7|R/2 = do, which, together with equations (5) and (13), enables us 





to write the equation for r as 


=~ 
7 


1 
f = tom-*leinthMy- coshtetanhe(2*]]? (41) 
, | F do 


io 


TPo 


where (42) 


~ 2HUB,,. 
Assuming ¢/¢) = 8/8), and putting o 1 (equation (38)), we find from 
equations (2), (41), and (42) that for the example of section 5.11 


»\1 
2 exp! ir tan sinh cosh? 1 tanha(*)?, (43) 
l | 8) } 
and = oH T\. 


When 7 has the value —7z/3, these equations yield the values given in 
columns 5 and 6, Table 2. A comparison of columns 4 and 6 shows that 
the maximum error of equation (41) for this example does not exceed 
6 per cent. of the velocity increment. 

The phenomena of jet adherence to curved bodies has been investigated 
experimentally by Squire (3). If a theory—perhaps a boundary-layer 
theory—could be established to predict 7, the angle of jet deflexion, then 
the method given above could be used to calculate the velocity distribution 
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over the body. It is worth pointing out, however, that if s is the distance 
measured along the surface from a point of separation, the inviscid theory 
yields in general a singularity in pressure gradient of the type s-} (cf. 
equation (43)). This leads to obvious difficulties in the application of 
boundary-layer theory to the calculation of separation points using the 
inviscid velocity distribution resulting from some assumed value of 7. 


5.2. The design of a bend in a channel 

As a second example consider the flow in the channel shown in Fig. 4 
To avoid boundary-layer separation, the inner wall A B is designed to be a 
constant pressure wall with boundary-layer suction applied at A. 


U7 








The w-plane is similar to the one shown in Fig. 2, and the solution can 
be obtained by replacing 4, by ¢)+7h in equation ( ma In this example 6, 
iszeroin —0o0 < d < O, igi? is equal tov in 0 < ¢ < &@, while 7, and r,, are 


constants. Thus i (34) yields 


2r cosh z w)/2h)\4 
f(w) = r,+i7 + — tanh 1} (Fo—w)/2h\$ (44) 
a \cosh 7(¢,+ w)/2h] 
a , 
where rs [i Bd(log - } 
q a 
The condition f = 0 at w x leads to the relation 
TT ‘ r bd 
tanh Po sech “4, (45) 
2h T 


which implies that given H, U’, andr, V (the velocity on A B) and the length 
of AB cannot be assigned independently. The shape of the free surface 
follows from equation (44) by putting w = ¢+7th, which results in 


7 _, 27 , {sinh 7(¢y—¢)/2h\ 3 
Ae) = 7 ee ee ee * 


which can be reduced with the aid of (45) to 


tanh(z¢/2h) sech(zr,/7)cos(76/7). 
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On the free surface 4 = Vs, and so using (5) we have finally the intrinsic 
equation for AB: 
tanh(7Vs/2H UB.) sech(z7,/7)cos(7@/7). (46) 
The total length of AB is plainly given by 
28s) = (4HUB,,/7V )tanh-'sech(zr,/7). 

5.3. The flow of a stream up a step, with boundary-layer separation 

As in Fig. 5, suppose that a stream flowing along a straight wall meets 
a step HC; then the boundary layer will separate at some point A in front 





o neeee= 
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of E. We shall assume that the separated streamline AC is one of constant 
pressure up to C, where there is a sudden drop in pressure to the value at 
infinity on CF. 

The w-plane for this example is as shown in Fig. 1 (c), and hence (see 
comments in section 2) the solution can be obtained directly from equation 
(34) by replacing ¢ by 6—d¢p, and letting ¢, tend to infinity. It will be found 
that 


0 
om ee 2 f ,. { 1—exp(—aw/h) \3 
f(w) r,+76,-4 -f | tanh Ni] aa. ~d*/h)| d0,(b*) 

a. - miei 

; | tan-1/ CXP(— 7 mt 12 dra(b*) 
T . \l1—exp(—7¢*/h)} 

d*=0 

at Ge {exp(—a7w/h)—1)3,. . 
= | tan ‘end nd* h) i| dr, (*). (47) 


In the example of Fig. 5, #, = 0, and if the pressure discontinuity at C, 
where ¢ = a, results in a discontinuity of o in 7», then (47) becomes 
») , 9 1 
: 20 exp(—7rw/h)—1\4 
f(w) = tan ( 1". 


|1—exp(—za/h) 


7 


At w = —o, f = 0, and thus o = 7, i.e. 
ail fy 34, (me mw/h)— 1\3) 48) 
Fw) A) = Ue exp(—za/h)] | 
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If s, is the length of the curve AC, then clearly 
a V 8. 
The lengths of AF and EC are given by the real and imaginary parts of 
. ; if ; 
I exp(2@) ds - y | exp(2@) dd, 
0 0 
respectively . On Vo od a, from (48) 
2 (3 exp(—7¢/h ) 
\) 


0 4Atanh — 49 
TT l—exp(—za/h (49) 


from which it follows that 


] |] —e-maih [ e’ tanh(76/2r,) dé 
Vr, cosh?(70/2r,)—(1—e-7@")sinh?(7 2r,) 


v0 
The imaginary part of this integral can be obtained in an algebraic form. 


If 1/2 is the length of EC, it is found that 


2h 2r 
l =e cosech ? { cosh| A tanh rl exp( Tra h))). (50) 


7 


An application to boundary-layer theory of the special cases of equations 
) 


(48) and (50) when h © and B, | (incompressible flow in an infinitely 


deep stream) was made in a recent note by Lighthill (2). 


5.4. Helmholtz flow of a jet impinging normally on a flat plate 
If the boundary layer of the previous example separates only at C, then 


the Helmholtz flow shown in Fig. 6 will result. 





The solution can be written down directly from equation (47)—there is 
a single discontinuity in 0, of 47 at ¢ —a, and hence 
= 1—exp(—azw/h)\} - 
/(w) lia +tanh-! 4 ie (51) 
” l—exp(za/h) 
It only remains to calculate the length of EC, which can be done analyti- 
cally in the case of incompressible flow. For this case it readily follows from 
equations (8), (10), and (51) that if 1/2 is the length of EC, and 


p exp(za/h)—1, 
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l hf by ogi i: +{p/(p+1)}) , 
-= ' 7 = “ti. Re 
aes 2 6h | (FI S\7- {p (p+) yyy al | (52) 


EC can be regarded as one half of a flat plate placed midway in an 
impinging jet. With the aid of Bernoulli’s theorem, we find that the drag 
on this plate is given by 





0 


D = pU? | testa dp 


( ( l —e-2"\er dd 
° 
= 2pU | sinhr dd. 


—a 


I 


The value of r on EC follows from equation (51); after a little algebra there 


results p= 2hpl "( ] —e-ma/2h) (53) 


For small values of the ratio //H, equations (52) and (53) yield 


+a pun) ( nm 1 
“\V44+r]\ 3(44+7)2 Ay’ 


the first term of which is the well-known classical result for a deep stream 
((7), section 12.21). 


5.5. A generalization of Borda’s mouthpiece 

In the ‘generalized Borda’s mouthpiece’ shown in Fig. 7 the flow comes 
from the region outside the tube and separates at A and B to form a free jet. 

The w-plane is of the type shown in Fig. 1 (d). Equation (14) yields for 
this case (¢, = 0) 
{ tanh(zw/2h) \4 


a ae momma chad heed J 
expt ‘) \tanh(zw*/2h)) ’ 





where w* is either 4* or 6*-+ ih depending on the range of integration. In 
this problem r, = 7, = 0; thus (27) can be written 


0 


, ei ae P tanh(zw/2h) \4 
7) 0 j= té } 1| iinindniatntmdinnias he i *\__ 
fw) wAT 7 | ei \tanh(7* 2h)) d80($*) 
o*=—o 
0 h h 

2 tanh(mw/2h) )4 . 

ute tanl _1{ tanh(zw/2h) \3 * 54 
7 | _ \coth(7*/2h) ; d0,(6*), (54) 


Pose 
where 6, and 6, are the values of @ on w = 0 and w = th respectively. 
One interesting feature distinguishes this example from those treated 








COMP 


above 
in Fig 


to f(u 
term 
show] 
descr 


me 
Di 


wh 


be 


here 


Cam 


mes 
jet. 


ior 


In 








or section 3. 





COMPRESSIBLE SUBSONIC FLOW IN TWO-DIMENSIONAL CHANNELS 281 


above. It is that there is a discontinuity of 27 in @ at 6 = —oo, ie. at D, 
in Fig. 7, which, from (54), contributes a term 

4 tanh-li(tanh(zw/2h))! 
to f(w). In the simple straight-sided Borda’s mouthpiece this is the only 
term in f(w). In the more general case, when @ varies continuously as 
shown in Fig. 7, the integrals in (54) must be evaluated by iteration as 
described in section 4. 








Doo @= 277 oe | 
B 
fle th aa ie at antl 5) 
01927 se <==}9-0 < 
=O * 
A 
D | 
a G=0 al 
Fic. 7 


For the straight-sided mouthpiecet 
f(w) = 4¢tan-)(tanh(aw/2h))!. (55) 
The intrinsic equation of the free surface, 4 = 0, 0 < ¢ < o, on which 
Us = 4, follows directly from (5) and (55). It is 
6 = 4tan-1(tanh(zs/2H))', (56) 


which is the classical result (7). 
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APPENDIX 
The exact equations (8) 
fale) oL 26 eL 
-—_(1— M*)— = 0, —+— = 0, 
on ag as * en 


where L — log(U/q), when transformed to the (¢, y)-plane defined by equation (1) 


becom«e 06 Poy v2 oL 0 26 p oL 
2 — a <= = Vs, al — op, 2 
op p ed Op © po OY 
1.€. from equations (2) and (3) 
o6 or 26 61 @r 
— — = 0, — —— es 0. (57 
ous ” Ob apt m dys ’) 


} Although f+ c as ¢—> ©, it is easily verified that f satisfies condition (a) of the theorem 
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Equation (3) can be written 


— 1/(y—-1) 
m = (142 at)” (me) 


4 


- as. M'+0(M*), 


where y is the ratio of the specific heats. Thus 





aa y+l 4 4\1 6 6 
Mm = Mo+ 5 (M*— M4,)+0(M*— M3‘), 
and so for subsonic flow it is a fair approximation to write m = ma. Then (57) can 
be written PY) Be 
¢ er 


- ae 06 —0 
Ama) Ob , eb * (mab) ‘ 


from which it follows that f = r+7@ is an analytic function of w = d+ima wp. 
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{ NOTE ON A PAPER BY DAVIES AND WALTERS 

ON ‘THE EFFECT OF FINITE WIDTH OF AREA ON 

THE RATE OF EVAPORATION INTO A TURBULENT 
ATMOSPHERE ’F 

By L. E. PAYNE (Institute for Fluid Dynamics and Applied Mathematics, 

University of Maryland) 

Received 27 November 1952] 


SUMMARY 


ie solution given by Davies and Walters for the concentration of water vapour 


| 
tany point in a turbulent atmosphere above an infinite rectangular saturated area 

y be transformed into a form which is more suitable for computational purposes. 
lhe integral expression for the solution given by the authors is replaced by a rapidly 


nverging hypergeometric series. 


Iv a recent paper (1) Davies and Walters were concerned with obtaining 
solution to the equation 
CO*y 0%) (l—m) ey 
x X 4 2 nt, ¢ >t (1) 
such that y vanishes at infinity and satisfies on the y-axis the conditions 
(a) X/ Xo 1, for ni SC) 


Du = (2) 
(h) X 0. n > € 


ippears from their discussion that the authors assume the inequality 
)<m <1. This assumption is made in view of the available experimental 
data. It should also be noted that their formula (4.5) would be meaning- 
ess without such a restriction. 

The equation (1) falls into the class of equations which were considered 
by Weinstein (2). In the terminology of Weinstein, y is an axially sym- 
metric potential defined in the meridian (yn, z) plane of a fictitious space of 

mdimensions.{ The particular solution to (1) which satisfies (2) is merely 
the electrostatic potential of a finite wire in the 3—m space. In order to 
solve this problem Davies and Walters introduced a system of prolate 
spheroidal coordinates and sought a solution to the differential equation 


O*y  d*yx Ox CX 
x X 1 (1—m)coth B X 4 (1—m)ecot 0°X — 0, (3) 
CD*~ c= C B C 6 
This wo vas performed under the sponsorship of the Office of Naval Research. 
The terminology is justified by the remark that if we replaced the factor (l1—m) by 1, 
vould be axially symmetric potential in a 3 space. This analogy to the 3 space has 
ved useful in various investigations. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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with the conditions that 


(a) x/xo>lasB>O ) 
(b) x>0 as B>o 
; (4) 
ox 
Cc “4 >50)atd=O0,7 
(°) 06 
The authors stated on page 473 that the Legendre P and Q functions which 
are the general solutions to (3) cannot be so chosen as to satisfy conditions 
(4). We shall show, however, that for 0 < m < 1 it is possible to express 
the solution to the boundary value problem in terms of the Legendre 
function Q—}” (see equation (7)). We shall show further that there exists no 
solution to the problem if (l1—m) > 1. 
A general solution to (3) which for m < 1 satisfies the second of con- 
ditions (4) (see Hobson (3), p. 234) may be expressed in terms of functions 
of the type 


sinh?” sin?”6 Q; *"(cosh B)[A,, P;, *”"(cos 6)+ B, Q; "(cos 6)], (5) 
where the P and Q functions are the Legendre associated functions of 
general type. Now for m < 1 and 0 < @ < z, as is the case in the problem 
under consideration, the following identity exists: 


cos 4mm P—}"(cos 0) + (2/7)Q—}%(cos 8) 


Sm 





2- 7 (1—m)(sin @)-3” = 
= - a sin ys)-”" d 6 
Terqd—my Jone 
0 

(see Hobson (3), p. 250). Thus it is apparent that a solution to (3) which 

satisfies sialieics 4) for 0 << m < 1 is given by 

2-Mm—2Zehmri 

sinh?” Q—}""(cosh B). (7) 


xiXe = Td—m)T(4 m) 


This expression can be shown to reduce to equation (4.5) in Davies and 
Walters’s paper. We merely make use of an integral representation for the 
Q function (see Hobson (3), p. 257) and write equation (7) as 


slog (cosh B 4 1)/(cosh B—1) 











(3) 

. sinh w)-™ dw 

X/Xo rs Ga— m)\ TQ ‘™m) } ( in w) aw 
0 g 
coth—(cosh B) \ 

21 (3) i , 
= aa | (sinh w)-” dw 
I{3(1—m)}T(4m) 


We now make the substitution sinh w = cosech¢ and note that 


cosh w = coth ¢. 
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Thus (7) reduces to the form 


21 (4) [ ‘ : 
a RS ‘Sa s I m—1 d : 9 
X/ Xo ri 1 m)\T (gm) ‘ (sin 1g) d, (9) 
B 


which is equivalent to Davies and Walters’s formula. Equation (7) may 
of course be expressed in series form as 
(4)cosh”-48 


1 
5 
Pf{4(3—m) tT (4m) 


F{4(2—m), 3(1—m); 4(3—m), sech?f]. 


X/ Xo 
(10) 


This form of the solution is more suitable for computational purposes than 
the integral representation given by Davies and Walters. 

It is obvious on physical grounds alone that the value of m in Davies 
and Walters’s paper must lie in the range 0 < m < 1. However, it is of 
mathematical interest to note that if 1—m > 1 there exists no solution 
to equation (1) which satisfies the prescribed boundary conditions. 

To prove this we observe that any solution y{1—m} of (1) may be ex- 
pressed in the following way (see Weinstein (2)) : 


x{l—m} = z™y{1+m}. (11) 


Equation (11) shows that to each function y{l+m} corresponds a 
uniquely defined (up to a multiplicative constant) function y{1—m} and 
vice versa. For m > 0 let us assume that there exists a solution to (1) with 
the factor (1—m) replaced by (1+m) such that conditions (2) and the 
condition at infinity are satisfied. Then according to (11) there exists 
a function y{l1—m} which vanishes on the line z = 0 and at infinity 
((22+- »?)'t+™yf1 +m} = const.). But a function y{1—m} which is analytic 
throughout the region exterior to the line 8 = 0 in the upper half-plane 
z >0, and satisfies equation (1) at every point in this region must attain 
its maximum or minimum on the boundary. Therefore y{1—m} and hence 
x{l+-m} must be identically zero. This means that no solution to the 
problem exists for (l1—m) > 1. The case m = 0 corresponds to the three- 
dimensional electrostatics problem of a finite wire for which it is known 
there is no solution. 

On the other hand, it is apparent that solutions to the problem do exist 
at least for some values of m in the range m > 1. We need merely make 
use of the relationship 


i 


a ’ 
oxi l—m 
xi, —1—m; nm-! <i 


— + x{l—m}, (12) 


on 


which permits us to obtain a solution to (1) when 2 < m < 3. 
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For m in this range the solution may be written as 


f Y 
xi 1—m} 





r'(4) } (sinh B)”"-? cosh B sin? es (m—2) (sinh 6)”-3 dé 
= = 35 —~— | (s 3 dd). 
I'{3(3—m)}T(4m) sinh?8-+ sin20 2 
B 
(13) 
Similar treatment may be employed to obtain solutions for all positive 
values of m in the ranges 2n < m < 2n+1. The remaining values of m, 
though perhaps of some mathematical importance, will not be considered 
in this paper. 
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TORSION OF A SEMI-INFINITE BODY AND A 
LARGE THICK PLATE 
By YI-YUAN YU (Washington University, St. Louis)+ 
[Received 6 January 1953] 


SUMMARY 

Michell derived the equation for torsion of bars of varying circular cross-section. 
Using this equation, and by means of Fcurier—Bessel series, Purser analysed the local 
effect near the ends of the bar due to an arbitrary distribution of torsional shearing 
forces over its terminal sections. A similar analysis is given in this paper for the 
torsional problem of a semi-infinite body and a large thick plate, which are essentially 
bars having very large radii. The Fourier—Bessel series is replaced by a Fourier— 
Bessel integral. 

Stress and displacement components for the general cases are expressed in integral 
forms. These are then used to solve particular problems. In the first problem the 
torsional shearing force is distributed linearly within a circular area at the boundary 
of a semi-infinite body. The second problem consists of a semi-infinite body sub- 
jected to a concentrated torque at its boundary, which is considered as the limiting 
case of uniform shearing forces distributed over a finite circular area of the boundary. 
The similar problem for a large thick plate subjected to two equal and opposite 
concentrated torques on its surfaces is also discussed. Numerical results are given. 


1. Introduction 

For bars of varying circular cross-section under the action of torsion, 
Michell (1) showed that the displacement of any point is directed at right 
angles to the axial plane passing through the point. If this displacement 
is denoted by v, and polar coordinates (r, #,z) are used with the axis of the 
bar as the z-axis, the only two non-vanishing stress components are the 


cv 6 ov 
70 = Hl: -‘), “ee. (1) 


cor Oz 


shearing stresses 


u being the shear modulus. These stress components can be expressed 

in terms of a stress-function ¢ by the equations 

pe od 
- ’ 


r* Oz 


“r 


Hop (2) 


Ta- - — » 
Os r2 dr 


and the stress function satisfies the following differential equation: 


Hh 306 OH 9 


(3) 


Or? ror’ @ 


Using equation (3), Purser (2) gave an analysis of ‘local perturbations’ 
Syracuse University, Syracuse. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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due to an arbitrary distribution of torsional shearing forces over the 
terminal sections of a bar of finite size. In this paper we discuss similar 
problems for a semi-infinite elastic body and a large thick plate for which 
torsional shearing forces are prescribed on their bounding surfaces. (A large 
thick plate here refers to one of finite thickness but otherwise unlimited in 
extent.) The applied shearing force is expressed in terms of a Fourier—Bessel 
integral instead of a series, the latter being used for the case of a bar of 
finite size. 

After this paper had been written, it was found that the results given here 
as equations (6), (7), and (8) had been derived earlier in Reissner’s paper 
on torsional vibration of a semi-infinite body without using the stress- 
function (3). These results (3) also involve the time as an independent 
variable. However, the connexion of the static problem with the familiar 
torsional equation (3), the calculation of the stress-component 7,9, which 
is hardly discussed in (3) but is rather important in stress analysis, as well 
as some numerical results to be given here, seem to be of general interest. 


2. Torsion of a semi-infinite body 

The body will be chosen to be bounded by the plane z = 0 with the 
positive z-direction pointing into the body. At the boundary of the body 
torsional shearing forces are prescribed with their axis coincident with the 
Z-axis. 

Since the stress components 7,9 and 79, and displacement v must vanish 
in the remote part of the body, we may assume a solution of equation (3) for 
the present problem to be of the form 


b(r,z) = | Reda. 


0 


The function R = R(r) must evidently satisfy the ordinary differential 


equation @R 3dR 
tiny 4 
dr? or dr — ; ” 
the solution of which is R = Ar*J,(Ar), 


in which the coefficient A is independent of r and J, is the Bessel function 
of the first kind and the second order. Substituting for R in the above 
integral for ¢ and then for ¢ in equations (2) yields 


7,9 = | wANJ,(Ar)e-* dr, (5) 


0 


t92 = | pA, (Ar)e-> dA. 


0 
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The displacement v can be obtained by integrating equations (1); the result 
gives - 
vy = — | Ad, (Ar)e-* da. (7) 
0 
The coefficient A is in general a function of the parameter A and may be 
letermined from the boundary condition 
(T92)2-0 Fir), 
Fir) being a given function of r. If we express this in the form of a Fourier— 
Bessel integral as P 
F(r) | AJ, (Ar) dA y F(y) J, (Ay) dy 
0 0 
and compare it with the expression 
9 = | #AdAAS, (Ar) da 
0 
3 derived from equation (6), we obtain immediately 


x 
. 


A | y Fly) (dy) dy. (8) 
b 


« 
0 


3. Linear distribution of torsional shearing forces within a circular 
area at the boundary of a semi-infinite body 
Consider a simple example in which the shearing stress tg, at the boundary 
cts within a circular area of radius a and is given by the function 
(Ar (O<r<a) 


F YT) 
10 (a<r), 


(9) 


K being a constant. Substituting this into equation (8) yields 


r 9 


@" J,(Aa). 


if 
A Ky?J,(Ay) d 
2 | LY 14 y) ay yA 


By equations (5)—(7) the stress and displacement components are therefore 
Trp Ka? | J,(Aa)J,(Ar)e cA dy, 
0 


rp. — Ka2 | J,(Aa)J,(Ar)e- dA, 


' Ka" {1 a(ra)d,(Arje- da. 
BJA” 


0 
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Integrals of the type occurring here may be evaluated by a method 
developed by Nagaoka and described in detail by Terazawa (4). However, 
here we shall consider only the values of the stresses and displacement at 
the boundary, where the preceding expressions reduce to 


(7,9)2-9 = Ka? [ J,(Aa)J,(Ar) dA, 


0 





(9z)2-9 = Ka? | J,(Aa)J,(Ar) dA, 
0 
(v).-9 = — ied | t y,(da)¥J, (Ar) dx. 
‘ we lA 
0 


The integrals in these formulae belong to the type of the discontinuous 
integral of Weber and Schafheitlin. The first integral may be evaluated 
easily by means of Gubler’s result (5), giving the following two expressions 


for (7,9)2-0' _. (r\? —f. y2 
2Ka() A(b 3:5) (r <a) 
a" *Y a? 


3 2 
\ 7 ta 


The usual notation ,/,(a,b;c;z) has been used to denote a hypergeometric 
function. The integral is discontinuous at r = a and the stress (7,9),_, there 
is infinite. The integral for (7,)..) may be evaluated according to a formula 
given by Watson (6) for a special case of the discontinuous integral; the 


toh 


(7,9)- ee 4 


ohn 


result can be identified with the boundary condition given in equation 
(9). Finally, the integral in the expression for (v),_, may be evaluated by 
means of a result due to Sonine and Schafheitlin (7); thus 


r 


| Ka? M 2 

| Jyh (2, -4:2:3) r= @) 
=p \@ a* 

| Ka®/a\? _ 2 

a (" 2 (2453:5) —— 

OM F} r 


Numerical values of (7,g)..) and (v),_) have been computed and plotted 
as shown in Figs. 1 and 2. All the above hypergeometric series converge 
rather rapidly so that the first few terms of each of the series are suft- 
cient for the purpose of calculating most of the stress and displacement 
values. It is interesting to note that while (7,9)--9 18 discontinuous at r = 4, 


(v), a 


(v),-9 is continuous everywhere, and that the maximum value of (v),_, does 
not occur at r = a but at some point within the loaded area. 
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4. A concentrated torque on a semi-infinite body 
dlotted The problem of a concentrated twisting torque acting at the boundary 
nverge of a semi-infinite body with the plane of the torque parallel to the latter 
» suf may be solved by considering it as the limiting case of one in which shearing 
ement lorces are uniformly distributed within a finite circular area of the boundary. 
r=a The function F(r) for the latter case takes the form 
9 does ; (K (0<r<a) 


Fir) 


{0 (a<r). 
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The resultant torque M is equal to 

M = | K(2zr)r dr = $Kza’, 

0 
, 3M 

and therefore K = ——., 
27a 
Substituting this value of K into F(r) and the result into equation (8), we 
find 





3M f 
4 = omna® | y J, (Ay) dy. 


0 
The integration may be carried out by expanding the Bessel function into 
a series of ascending powers of (Ay) and integrating the result term by 


term; thus eS 
3M Cc “ —])r /Jay\2n41 
27a? | n!(n+1)!\ 2 
0 


n=0 





2 2n+3 


3M = (— 1)” A\2n+1  @2n 
27 n!(n+1)! 


n=0 
The limiting case in which the radius a of the circle approaches zero gives 
the case of a concentrated torque of magnitude M acting at the origin, for 
which the value of A is 


The stress components are, by virtue of equations (5) and (6), 


T= | a r2J,(Ar)e-*A dA = ee , (10) 
Jj * 4n(r?+-z?)! 
0 
r M.. 3Mrz 

T9, = [ 4, Mare <A d\ = meta’ (11) 


0 


the integrals having been evaluated by means of Gegenbauer’s formulae (8). 
The displacement v is similarly found to be 


SAO — (12) 
4a 4r(r?+-2?)! 


0 








It can be verified easily that the boundary condition is satisfied, for 79, 
vanishes at the surface everywhere except at the origin, where it becomes 
infinitely large because of the action of the concentrated torque. 
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OF 
= 
re 6+ F779) z- /M 
+ | 
4-L 
t rélz- of M 
) A 1 
= 0°5 1:0 15 20 
s=7/c 
Fic. 3. 
sf 
S 4.5L JAE OW 2 gf/M 
1 ~ 47 p10(V).--0/M 
5b 
( - t 
)-5 1-0 15 20 
$=7/e 
Fic. 4. 


Introducing a length parameter c and s = r/c, we obtain, by equations 


10) and (12), the values of 7, and v at the boundary, 


4rc¥(7,9):-0 _ 3 


i 3 

’ * | (13) 
4z.c*(v). o 1) 

M gt 


lhe dimensionless quantities on the left-hand sides of these equations are 
computed and plotted against s as shown in Figs. 3 and 4, in which similar 
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curves for a plate problem to be considered later are also plotted for 
comparison. 


5. Torsion of a large thick plate 

The foregoing analysis for a semi-infinite body may also be extended to 
the case of a large thick plate with prescribed torsional shearing forces 
acting on the two bounding planes z = +c, where the half-thickness ¢ of 
the plate can have any finite value. The boundary conditions are 


(T9:)- 1 = Fir), 
(T9:):=-e — G(r), 
in which F and G are given functions of r. 
Without imposing the limitation that the stress and displacement com- 
ponents should diminish with increasing values of z, the solution of equation 
(3) obtained in section 2 must now be modified to take the following form: 


_ 


y= J (Re-A+- R’e*) dx, 


in which R and R’ are both Scliteeis of r. As both R and R’ must satisfy 
equation (4), we have 
R(r) = Ar?ZJ,(Ar), R'(r) = Br*J, (Ar), 
where B as well as A are constants independent of r. Therefore, 
¢ = - | ( (Ae-A+ Be*)r2J,(Ar) da, 
0 
from which we obtain 
0 = | p(Ae-— Be*)AJ,(Ar) da | 
j | 
to, = | p(Ae+ Be)AJ,(Ar) dd $. (14) 


0 


y= (—Ae-*A+ Be*)J, (Ar) dx 
0 J 
To determine the coefficients A and B, we express, as before, the given 
functions F and G@ in the form of Fourier—Bessel integrals: 


F(r) = i AJ, (Ar) da [ y F(y)J, (Ay) dy, 
0 0 


G(r) = | AJ,(Ar) dd | y Gy), (Ay) dy. 


0 0 
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The second of equations (14) then yields the following integrals for the 


houndary stresses: 


= 
(Tiilen te | (Ae cA + Berd) J, (Ar) da, 

0 
(raede--e = { p(Ae+ Be~')AJ,(Ar) dA, 


which are equal to F(r) and G(r) respectively. By a simple comparison 
between these two sets of expressions we obtain 


p(Ae-+ Be) ( y F(y)J,(Ay) dy 
: +, (15) 


p(Aer+ Be) = | y G(y)J,(Ay) dy 
0 
from which A and B may be obtained. 
For the symmetrical case in which F(r) and G(r) are identical, the stress- 
function ¢ must be an even function of z; therefore, we must have 


A=B. 
Equations 14) now become 
7,9 = — | 2uAAsinh2dJ,(Ar) dA 
0 
7 { 2uAAcosh2AJ,(Ar) dA, (16) 


"2A sinh 2A J,(Ar) dA 


and the two equations (15) reduce to one, namely 
2uA coshcdA = [| y F(y)J,(Ay) dy, (17) 
0 
from which A may be determined. 

Noting that (v)._, is zero by equation (16), we conclude that equations 
(16) and (17) also constitute the solution for a large plate of thickness c 
with one side subject to torsional forces and the other side fixed to a rigid 
base. 


6. Concentrated torques on the opposite bounding surfaces of a 
thick plate 
We shall consider the symmetrical case of two equal and opposite con- 
centrated torques of magnitude M acting on the bounding surfaces z = +c 
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of a thick plate, the line joining the two points of application passing through 
the origin. Following along similar lines as for the case for a semi-infinite 
body, we have, by equation (17), 


I 3M [ yn 
P a = M4 J, r ] pe 
. 2 cosh ca roast Ina? | yd, (Ay) dy = 
0 


87 cosh cA’ 


Substituting this value of A into equations (16), we obtain 





* M22 sinh zd 
Tré — — 


4m cosh cA 


J,(Ar) da, (18) 
0 
* MA? cosh zA 


iil . 47 coshcA al ae 


* Mixsinh 2A 
47 cosh cA 


0 


J, (Ar) da. (19) 


We shall investigate the values of these stress and displacement components 
at one of the two bounding surfaces of the plate, for instance at z = —c. 
where the above expressions reduce to the following: 


4)..-0 = = A* tanh cd J,(Ar) da, (20) 
. 0 
Mf, 
(Ton)eu-e = =I AJ, (Ar) dA, (21) 
i 0 
M 
(v),--. = { Atanh cAJ,(Ar) dA. (22) 
doe 


0 
That (79,),--, is zero may be seen by comparing its integral expression 
in equation (21) with the integral in equation (11); thus, using the latter, 
we have 
- Bre 
| A*,(Ar) dA = lim | A*J,(Ar)e~ dA = lim —""* __ — 0, 
4 0 9 230 (r?-+2?)3 
which is true everywhere except at » = 0, where the concentrated torque 
is applied and where the stress becomes infinite. 
Values of (7,g).-_, and (v),__, may be obtained by numerical integration. 
We first make the substitutions 


cA=2z and s= r/c, 
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so that equations (20) and (22) become 


Mf, 
(7,9). 3 | x* tanh xJ,(sx) dx, 
° 0 
Mf 
(v),--, - | «tanh x J,(sx) dx. 
: 4zrpc* | 
, 0 
Noting that, for 2 > 5, 
tanh x 1-0000 to four decimal places, 


we may write the above expressions with sufficient accuracy as 





= a ? x* tanh xJ,(sx) dx — | x*J,(sx) dx + | x*J,(sx) dx, 
0 0 0 
as ‘ Jeane «tanh x J,(sx) dx — | xJd,(sx)dx+ | xJ,(sx) dx. 
0 0 0 


The following relations are obtained if we put z equal to zero in equations 
(10) and (12): 


x 


€ 


“7 3 
| 2J,(sx) dx = 


a? 
3 





0 
° 1 
| xJ,(sx) dx = 3 


Therefore. a 


47c*(7,9) , - } 
— —=+ «*(tanh x—1)J,(sx) dx | 
4 s | 
(i a (23) 
dorc?(v ] r 
BOW) e= = + a(tanh x—1)J,(sa) dx 
M ee 


Using Simpson’s rule and an interval of x = 4, the numerical values of the 
definite integrals in these equations can be computed for different assigned 
values of s. The final results are tpbulated below: 





1/s? 470(T19)2---/M 4mpc*(v),-_,/M 





x x a 
2°7778 | 13°8374 2°6588 
I*O “0000 1°0000 | 2°8974 0°8416 
o*5102 0°9g042 0°3433 
c 75 0*2500 02550 o*1065 





Values in the last two columns plotted against the corresponding values 
of s are shown in Figs. 3 and 4. The resulting curves are lower than those 
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for the similar problem of a semi-infinite body considered earlier. Com- 
paring equations (13) with (23), we see that the difference is due to the 
additional integral terms in the latter.t The difference may be seen by 
comparing the second and third columns respectively with the fourth and 
fifth in the above table, values in the former being also equal to the 
dimensionless ratios on the left-hand sides of equations (13). 
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t+ To be exact, these integrals should have an infinite upper limit. 











AN 


By 


in i 
por 


of 1 


(0 the 
en by 
h and 


O the 


‘ditior 


lrch 


under 


(Mac. 





Com- 





—— 
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SUMMARY 
A variational method is developed for the calculation of the propagation constant 
in inhomogeneously filled wave-guides which is applicable to waves whose field com- 
ponents are given in terms of one scalar potential. The conditions for the existence 
of pure TE or pure TM waves in such a wave-guide are also considered. 


1. Introduction 
\ECENT developments in electronics, such as the linear accelerator, have 
necessitated the use of wave-guides which are partly filled with dielectric. 
As the object of inserting dielectric into wave-guides is to reduce the phase 
velocity of waves in the guide, it is necessary to find the propagation con- 
stants of waves in such a guide. The calculation of propagation constants 
associated with a wave-guide which is filled with material such that the 
electrical constants vary across the guide, but not along it, involves in effect 
the matching of fields which satisfy Maxwell’s equations in the various 
media at the interfaces between them. Even when the system is simple 
geometrically, this involves heavy numerical work arising out of the solution 
ofa transcendental equation. In order to avoid this, recourse must be made 
to approximate methods. In general the modes which can exist on such 
structures are of a complex character and it is not usually possible to reduce 
them either to modes without a longitudinal electric field component (trans- 
verse electric or TE type) or to modes without a longitudinal magnetic field 
component (transverse magnetic or TM type). 

A full discussion of all the work which has been carried out on this 
subject has been published elsewhere (1). 


2. General properties of propagation on cylindrical structures 

In order to discuss such a system it is convenient to consider first the 
electromagnetic fields which can exist in an isotropic homogeneous loss-free 
medium of dielectric constant « and permeability «. The Giorgi system of 
units will be used throughout the course of this paper. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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Assume that the walls of the wave-guide are parallel to the z-axis and 
let u,, U, together with z form a right-handed orthogonal system of co- 
ordinates; then an infinitesimal line element is given by 


dr = i,h, du, + i,hy dug + i, dz, 


where i,, i,, i, are unit vectors in the direction of u,, u., z increasing. 

It may be shown (2) that in such a medium it is possible to express any 
electromagnetic field in terms of two Hertzian potentials Ili, and 1%, 
where II, I1* satisfy the wave equation. (II generates a pure TM wave and 
II* generates a pure TE wave as indicated under.) 

If now there are several contiguous cylinders of such media, a necessary 
and sufficient set of conditions for the matching of the fields at the interfaces 
is the continuity of L,, E., H,, and H., where E, and H, are the components 
of electric and magnetic fields respectively which are perpendicular to the 
z-axis and tangential to the surface of discontinuity. Thus, in general, the 
field is defined in terms of two potentials which are implicitly related through 
the continuity conditions at the interface. Under certain conditions, how- 
ever, it is possible to define the field in terms of one potential only, thereby 
simplifying matters considerably. 

Although in any physical system the value of the permeability does not 
vary appreciably from its value for a vacuum, it does not involve any severe 
complication in the analysis to consider varying as well as «, and accord- 
ingly this variation will be allowed for. 

In the event of there being a pure TE wave, the electromagnetic field in 
a homogeneous isotropic medium is given by 


, all * “2 aT s 
Kae, H, = B orl* 
hy Uy h, eu, 
gp eg FT i 
h, eu, “hy 0Uy 
E, = 0, H, = (k?—g2)I1*, 


where a factor exp{i(8z—wt)} is omitted and k? = wue. 

Suppose now that the electrical constants are functions of u, only. The 
equations of the interfaces are of the form wu, == C and the tangential fields 
E,and H, are now £, and H,. For a solution the following quantities must 
be continuous at the interfaces: 

E Lou ell* ip oll* 


es - Te fl 
2 


h, eu, “hy Oty 2) 


H, = (J®@—82)11*. 
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and } In effect the three quantities 
f CO- 
[I * oll* = m 
Pipe ——, (k®—p?)I1* 
| Uy OUs 


must be continuous at the interfaces. 
[I* satisfies the wave equation and this is of the second order so that two 


any boundary conditions are necessary and sufficient to determine I1*. Unless 














[*is, }  @fl*/éu, = 0 the problem is overdetermined and there is not a solution. 
and If, however, this condition is satisfied identically, sufficient conditions are 
given by the continuity of 
sary , , 
ae ye ad I. 
ents ss 
) the In an exactly similar manner, it may be shown that if the electromagnetic 
the constants are functions of uw, only, it is only possible for a TM mode to 
ugh propagate if cll /éu, 0. 
low- A further condition is imposed on TM modes through £, having to vanish 
‘eby on the wave-guide walls. As éIl/éu, 0 it follows that IT (and hence £,) 
does not vary along the interfaces (which are of the form u, = const.). 
not Thus if there is a metal boundary given by u, = constant, Z, must vanish 
vere everywhere and the TM modes cannot exist. Thus TM modes may exist in 
ord- guides such as a circular wave-guide, but cannot exist in guides such as the 
rectangular guide or the septate guide. 
din The analysis above is applicable to any cylindrical structure and it is not 
merely true for wave-guides. It follows therefore that on a cylindrical 
dielectric structure, waves which are transverse electric or transverse 
magnetic with respect to the axis of the cylinder can be propagated if there 
isno variation in field over the interfaces and if the system is ‘geometrically 
1) similar’. That is, the equations of the interfaces are of the form u = u, and 
the wave-guide (or cylindrical structure) occupies 
A Uy * Bp, U< Us « D. 
The 3. Energy equation for a guide with arbitrary dielectric filling 
elds The only case of practical interest (and indeed the only tractable one) is 
lust when everything is independent of z. Let a typical cross-section of the 


wave-guide be as indicated in Fig. 1. The cross-section is divided up into 
several regions such as p with electromagnetic constants y.,, €,,. 

Let t be a unit vector along any curve in the cross-section and let n be 
aunit vector along the normal to such a curve. Then the conditions at the 
walls of the guide, assumed perfectly conducting, are that E,, E,, H,, vanish. 
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The conditions at the interface between two media are that E,, H,, E., H 
are continuous and so (E/H).n is continuous across an interface and 
vanishes on a wall. 


It has been shown (3) that, for an isotropic homogeneous loss-free medium 

V.S* = hio(vwH.H—<E.E), 

where S* = 4E H (the complex Poynting vector) and the tilde indicates 
that the complex conjugate value is to be taken. 

Let z) be an arbitrary section of the guide and consider the integral of 

V .S* over the volume V, of the cylinder having the region p as cross-section 

and extending over z) < z < 2)+2z7/8. Then J, is a region containing a 


homogeneous isotropic medium. By Green’s theorem (4) we have 





where dV is an elementary volume and dA is an elementary area on the 


surface A, surrounding the volume JV, van be regarded as made up of | 
the following surfaces: 

(i) A,o—the portions of the planes 2) + 27/B intercepted by 
—the portion of A, which is in contact with the wave-guide 
A series of areas A is that part of the surface of region p 

in contact with region q. 
It follows that 


‘WH. Al—<E.B} dv 


S*.dA+ | S*.dA4 
yan 


[ S*.dA, and similarly for 
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Because 27/8 is a wavelength, the value of S* at z)+27/B will be the 


same as at z = z». The direction of dA at z = z,+27/B will, however, be 
exactly opposite to that of the corresponding dA at z = z. Thus I), 
Furthermore (E ( H).n vanishes on the wave-guide wall ond thus I. =z @. 
Hence, integrating over the whole cross-section of the guide and over 
<< 2% +2n/B, 
| (wWH.A—cE.E) dv > | (WH. A—<E.B) av SS hs 
J Vy iw p qa 
: . pF 
also lon S*.dA : —_ S*.dA —I,» 
A, in 


because for two contiguous regions the outward drawn unit normals are 
oppositely directed. Thus 


( (uH.A—<E.B} dv = 0, 


or | dz { {uH.A—cE.B} dw = 0, 

Ww 
where W is the cross-section of the wave-guide. Because all the field com- 
ponents are of the form F expf{i(wt—fz)}, the quantity (uwH.A—<«E.E) is 
independent of z and so the energy equation assumes the form 


| {ujJH|?—e|E|\*} dW = 0. (4) 
W 
Equation (4) expresses the fact that the mean magnetic energy is equal 
to the mean electrical energy. This is analogous to the result in mechanics 
where the mean kinetic energy in a vibrating system is equal to the mean 
potential energy. From this one is led to speculate as to whether it is 
possible to calculate propagation constants for wave-guides by approximate 
methods analogous to those used in determining the eigenfrequencies of a 
mechanical system. 
It will be observed that equation (4) has been derived under the assump- 


tion that £ is real, i.e. that the mode is propagating. The equation holds, 
however, equally well when the mode is evanescent. In this case f is 
imaginary and the integration takes place over an imaginary length of the 
guide, the result being formally the same. 

Now the electromagnetic field components are derived from two poten- 
tials, namely Il and II*. and, in general, the relation between II and II* 
will be complicated. Thus it will not usually be possible to express equation 
4) as a condition on one quantity only. It will be shown later in this paper 
that, if the field components are expressible in terms of one function only, 
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it is possible to make use of approximate methods, the quantity II (or I1*) 
occupying a position analogous to that of displacement in a mechanical 
system. 


4. Calculation of propagation constant 
A. TE waves in a geometrically similar structure 
Consider the geometrically similar structure composed of n materials 


such that the electromagnetic constants over a,_, <u, < a, are e,, p, 


(a) = A, a, = B). The other boundaries are given by u, = C, u, = D, 


m 


where C and D are constants. Then the TE field existing in the segment 
yy < Uy <a,, C < uy < Dis of the form 


Pp 
EH, = 0, Lune. 
h, Ou, 
I ee se ll H, = 0, (5) 
y h, ou, : 
E, = 0, H, = (ki—p*)I1. 


Thus the only non-vanishing electric field component is E, and it will be 
convenient to express everything in terms of this rather than IT (which 
obeys rather awkward continuity equations at the interfaces between the 
various media). 

Let us now assume that » and « are continuous variable functions of u, 
only. (The case of a discontinuous distribution can be accounted for by 
smoothing out the discontinuities in the usual manner.) 

Maxwell’s equations for any system are given by (5) 


V A E+ = oe 0, 
ot 
(6) 
VA = = 0 
ot 
Eliminating H we have 
VA (u1V A E)—w*cE = 0. 
Let E = fi, exp{i8z}. Then 
ey 1 @ e ‘ 
a pen h. ae k2 — 2b/u = 0. 7) 
Su, ah, he ou, 2”) ( Bb u \ 


The problem has thus been reduced to the solution of this equation subject 
to certain boundary conditions (4% = 0 when u, = A, B). 
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This is an ordinary Sturm—Liouville (6) problem, and an expression for 
# may be obtained by multiplying through by h,% and integrating over 
4<u,< B. It follows that 


du, 


. 
* Tkh, wv ] { ¢ 

: —— (kh, 
| | ws phy h. lou 2 wh 





B 


’ Re is 
| ~2 2 du, 
d F 


B. TM waves in a geometrically similar structure 
In exactly the same way as in section 4A, it follows from Maxwell’s 
equations, on eliminating E, that 
VA (e7?V A H)—ow*nuH = 0, 
and if H = wi,exp{i(8z—wt)}, it follows that 











| k2— 2 
. (how)| 4.2? .W 0. (9) 
6u,\eh, h. Gu, ~"'J € 
Multiplying through by h, % and integrating over A < u, < B the relation 
obtained for 8? is 
B 
i Lh, us l | 
a h. du,+P 
| | ‘ ch, h, {aur he | L 
A? — 4 - . (10) 
| ha du, 
P € 
4 
l @ silinal 
where h, (hs . 
' = | x hy eu, al. sic 
: da sti ] 
Now iwE = VAH — iBypi, + - (hg tb)iz 
h, rhs Cu, 
+7%4=B 
and so P hab Ey = 0 
ew] _, 
because EZ, vanishes on the walls u, = A, B. Hence 
» 
* [Rho l c 2 
| i i. (he) du, 
€ €, hy hg\eu, i | } 
Bp? 7 : (11) 
“he? 
| te du, 
‘i € 
A 


x 
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C. Propagation in a rectangular guide 
It has been shown that an inhomogeneously filled wave-guide cannot 
sustain either a pure TE or a pure TM mode unless that field is such that 
there is no variation in field along the interfaces (7). However, for a 
rectangular wave-guide such that the electrical constants only vary along 
one direction (Fig. 2 represents a typical guide 

















i of this nature), it is possible to split the field 
into modes such that either £,, vanishes or H. 

by a C20 pa vanishes. These are termed longitudinal seec- 
| tion electric (LSE) and longitudinal section 
ye ew magnetic (LSM)+ modes respectively. The calcu- 
Fic. 2 lation of the propagation constants for these 


modes depends as usual on the solution of a 
transcendental equation arising from the matching of the electromagnetic 
field at the interface. 

For simplicity only the case of an LSM wave (H, = 0) will be considered 
here. The result for an LSE wave can be derived in the same manner. 

Consider a rectangular wave-guide whose cross-section is 0 < x <a, 
0 < y < band whose walls are parallel to the z-axis and suppose that over 


Ly. < X < 4, (region p) the electromagnetic constants are .,,, €,, (% = 4, 
x, = a). Fig. 2 represents a typical guide of this nature. In a homogeneous 


isotropic region a possible field is given by 








ee: ae a 
cz" cy” 
, e e 
es mien (12) 
oxoy otoz 
os yp a aor 
; OXCZ : toy’ 
a2 a2 2 p2 
where u satisfies vf Pe in Pie zz @. 
Ox* oy o2z* ot? 


Now, in order that the boundary conditions on y = 0, y = b may be 


satisfied, 4 must be of the form 
= % ¢, sin. 
ig is Yn ‘ b 


n=1 





Consider, therefore, the field generated by 


. ny ; a 
Y = x(x)sin ™ exp{i(Bz—wt)}. (13) 

) 
+ Some confusion has arisen in (7) over the LSE, LSM notation. The wave with E, = ! 
has the E vector in the plane of longitudinal section and is thus LSE. This is consistent 
with the TE, TM notation. I am indebted to a referee for these remarks. 
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It follows, in order that % should satisfy the wave equation, that x must 


(eB |x = 0 (14) 
b? 


satisly dty 
dx? 


and the field components are given by 





;' -, ar . nm 
E, p 4 ; }x sin J H, 0, 
bh |’ b 4 

na dy niry etek . nrry 
E., - ONE ae. H —iwe iby sin —, 

Te” Sl uv pany 
ay . ni . nT na 
E ip<X sin “74 H, = iwe — x cos ~ J 

ax b a b J 


the exponential factor having been dropped. In particular a solution of this 
form is valid in region p provided that k*, yp, « are given their appropriate 


values. 


Multiplying equation (14) by e,, x and integrating over z,_,< «4% < a, 
. d*y (13 g2 n*7\ , Ia Q 
+ | k* —p?-— ——}y?| dx = 0. 
| “p| X da* r b? x 
Tp—) 
Integrating the first term by parts 
. > on Ip—O0 
‘dy\ P > nn \ , dy|” 
€, | x) 4 ia p? — —— )x® dx -le,Xx x 
| dx ; b? dx 
4 rr 


where the value at «,,-0 means the limiting value when x > x, from above 
or below respectively. It is not difficult to see that 


dy ° . . 7 
x2 & of the form GH, E,,, 


where G is a constant. 


It follows, summing over all p, that 


(15) 


0° 


Y [« | ~(3) +(8—e— | dx — G3 [H.E,y?~° 
a, ee Udy 14 
Since H., E,, are both continuous at x,, the product H, E, is also continuous 
there. In addition E, = 0 at x = a) = Oand « = 2,, = a. It follows that 
the right-hand side of equation (15) vanishes, and using the appropriate 
values of ¢, x, and k®, equation (15) may be written as 


, Us g2_ 3) _ (ZY | dx = 0. (16) 
L 
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If there is a continuous distribution, this is equivalent to an infinite 
number of regions p and the limit may be taken. Equation (16) is thus true 
also for a medium with continuously varying electromagnetic properties, 

It follows that 


- [ e[k?y?—(dy/dx)*] dx 
‘ p47 — 0 ' (17) 
' a 

| ex” dx 


0 





The case of LSM waves is formally identical and leads to 


a 


| [k2x2—(dy/dx)?] dx 











2.52 : 
B+—5 = 0 - : (18) 
| wx? dx 
0 
where H, = (B+. Joos “= x(aexpfi(B2—at)}. 


5. Expression of propagation constant in variational form 
It will be observed that in all cases (TE, TM, LSE, LSM) the propagation 
constant is given by an expression of the form 
P= 4-5 (19) 
Oi 
where A, B, C are positive definite (in other words, they vanish if and only 
if the potential (% or y) is identically zero). 

In order to discuss the matter further, it will be convenient to arrange 
that the quantity under consideration is always positive. The analysis 
which follows will be confined to LSM waves. The extension to other types 
of waves is very similar. 





9 9 9 9 9 9 na P 
Let x2 — k2.—k? and »? = ,,—p*~—- oe (20) 
where k?,,, is the maximum value of k°. 
It follows immediately that 
| e{«?x?+ (dy/da)} da 
oe = Q(x) (21) 
a R(x) 


[ ex? dx 


0o 


where Q(x), R(x) are positive definite. It is clear that y? > 0. 
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Consider now the variation of the equation 


(y?—k?)y? — | dz = 0, 


it being assumed that there is a continuous variation in dielectric constant 
and permeability. A discontinuous distribution is covered by regarding it 


* 


| < 


ny 


0 





sa rapidly varying continuous distribution. 
The condition that y* shall be insensitive to small variations dy in y when 
y+8y still satisfies the boundary conditions is given by 
a a 


P = dyd 
e(y?—k*)ydx dx — €- x! (dx) dx = 0. 
; dx dx 
0 0 
Integrating by parts and using the fact that 5y satisfies the boundary con- 
ditions this gives 


a 
“[d/[_ dx ae 
(< *| — €(y?—K?)x dy dx = 0. 
a dx dx 
0 
Now 6 x is, apart from its satisfying the boundary conditions, arbitrary. 
It follows immediately that y is that solution of the equation 





eae ox) + €(y?—x?)y = 0, (22) 
| which obeys the boundary conditions on the walls of the guide. 

The problem has thus been reduced to the Sturm—Liouville eigenvalue 
problem (6), and the following statements are true: 





A. There is an infinite set of eigenvalues y2 such that 0 < y? < y2 < ....+ 
B. There exists a complete normalized set of eigenfunctions {y,,(x)} such 
that 
] 
; al dx, 2 -2 ° 
i) é 1 e(y2—K?)y, = 0; 
dx\ dx | , 

li) x, satisfies the boundary conditions ; 
ili) } €y ro dx 0», (23) 


, 
: a X dy nS 
, m . ae ° " 9 
ai e| -K Xn Xm dx = Ym¥n Omns (24) 
\ da dx 
. 
0 
where 6,,,,, l if m = n and 6,,, = Oif m ~n. 
* Under certain cireumstances, two eigenvalues may become equal through degeneracy. 
example in an empty guide of rectangular cross-section with sides a and b, the eigen- 
alues are effectively f m?*/a*+-n?/b?, where m and n are integers. If a and b are com- 


mensurate, different sets (m,n) can give rise to the same f and the system is degenerate. 
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Thus, if G is any function which is piecewise continuous in 0 < x <q. 
it may be expanded in the form 


x 


G(x) > G,, Xn(*), 
1 


where G. = eGy,, dx. 








Let Ye = 


a 
| G? dx 


0 


then using the relations (23) and (24) it follows that 


xr 

2 92 
Gv 

2 n: a 


Y@ > 


tr 


1 

x 
¥G 
— 
n=1 


a 


Thus if G is an approximation to y,,, y?, is an approximation to y?, and 
an error in the approximating function generates an error of the second 
order only in y?. 


The discussion has up till now been in terms of the quantity 


yA 


which is always positive. It follows that B? (= k?,,,—y?—n??/b?) may be 
either positive or negative. The sequence y2, tends to infinity with m and 
thus there exist only a finite number of modes for which 8? > 0. That is 
to say there are only a finite number of modes which will propagate, all 
the others being evanescent. 

In general a wave-guide is used for transmitting signals and so it is 
necessary that only one mode should exist. That is, the frequency of the 
signal is above the critical frequency for the first mode and below the critical 


frequencies for all the other modes. These other modes are therefore | 


evanescent (see section 7). Thus it is only f? (or y?) which is of interest. 
? 


If now G is an approximation to the first mode, y?, will be an approximation 


to y} (in which the error is of the second order) and 


x x 

2 2 12/,2__ 2 
2 Gn ce > GilYn—Yi) 
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whence . P , n2a2 
Py “max ¥1— b2 
” > nn ak 
“max YG — he (26) 


Thus if G@ is an approximation to x, the value of A} is greater than the 
approximate value obtained and the difference is of the second order. 
It follows also that if several approximate functions G are taken, the best 
:pproximation is that which gives the greatest approximate value to fj. 


6. Application of Rayleigh Ritz methods 
Let y? be any one of the eigenvalues of 


d (_dx\ 


dx\ dzx| 


e(y* K*)x 0 (22) 
when y is subjected to the appropriate boundary conditions. Then it has 


been proved that 


| e{x*y?+ (dy/dx)*} dx 
, 2 0 es — = * (21) 


a 
9 
| ex? da 


Let ¢,(x) be some set of functions which is complete and orthogonal over 
<< a#< a and which satisfies the same boundary conditions as the y,,. 


Then y may be written in the form y = > 6,¢,. (A convenient set is that 
=1 


r 


composed of the modes which exist in the empty guide.) Consider the 
N 

approximation to N terms of x, x(N) > 6,¢,. Substituting x(V) in 
r=1 


equation (21) it follows that an approximation to y* is obtained of the 


form 


where the quadratic forms are clearly positive definite. 
Thus, y? is stationary for small variations, and if 
\ N 
F=) > (Q,,—y'?R,,)b,6, = 0 (28) 
1 s=1 


F /éb, ) for all s. That is to say 
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N 
and > (Q,s—y*R,.)b, = 9. (29) 
r=1 


The condition for this to have a non-trivial solution for all the 6, is that 
the determinant |Q,,—y?R,,| shall vanish. This gives rise to an equation 
of the Nth order in y? and so it is possible to find approximations for the 
first V propagation constants. 

It will be observed that the problem is formally identical with that of 
Rayleigh’s Principle in vibration theory and any discussion on such matters 
as the convergence of the y? to their appropriate values as N tends to infinity 
will be the same as in the mechanical case (8). 

It is possible to use relaxation methods (9) to obtain upper and lower 
limits to y* (and hence to £?) in exactly the same way as for the eigen- 
frequencies of a vibrating mechanical system. The energy equation for 
such a system can be written in the form 
N N 
Zz > Vier Us Vv 


eit sx ot i nt — —, (30) 


NN 
> 7 T 4, 
r=] s=] 
where the numerator and denominator are positive definite forms. The 
problem of obtaining bounds for y? is formally identical with that of obtain- 
ing bounds for w*. This can be done in the following manner. 
If G, = y?R,—Q, where 0R/eb, = R,, Q/eb, = Q,, then G, is analogous 
to the force F, corresponding to a displacement a,, Where 
ooT av 
=— q@-* a shoal aa 
éa, a, 


F, 


The analysis may be carried over by means of the following table. 


Mechanical system . ee T w? a, F, 
Wave-guide ‘ ‘ : 3 R y* b, G, 





N 
G Use 2 >, &, ' 
%+—(3) <"n<%— , (31) 
L > 5, R, 


where (G/R), denotes the algebraically largest value of (G,/R,) which 
corresponds with a positive product 6, R,, and y*, is the trial value for vi 
obtained from the relation 

2 = Q(x) 


R(x4) 


Whilst the interest has so far been centred on the fundamental mode, it 
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is possible by ‘subtracting’ this mode to deal with the second mode, and so 
on. This again is formally the same as in the calculation of the eigen- 
frequencies of a mechanical system. 


7. Critical frequencies 

In a homogeneously filled guide there corresponds to each mode a critical 
frequency below which it will not propagate (10). For this frequency 
92 — (. The same is true of inhomogeneously filled wave-guides. 

It will again be convenient to use LSM waves for the sake of illustration. 
In this case the problem is to find the eigenvalues of the equation 


d | dy) . . n2 772 

«—“\4 e[ k2—p2~- = 0 33 
dx\* dx| — )x is 
when 82 = 0 and y obeys the appropriate boundary conditions. In exactly 


the same way as for finding the propagation constant, it may be deduced 

that a 
( e[ (dy /da)? + (n*2?/b?) x] dx 

w? — oO ; (34) 


a 





. 212 d 
| pe xX” ax 


0 


which is of the form which has been discussed earlier in this paper. 


8. Application of method to wave-guides of arbitrary cross- 
section 
Consider a wave-guide of arbitrary cross-section filled with material of 
uniform electromagnetic properties and let % represent one of the rectangu- 
lar field components (for the present ys will not be defined further). Let x, y 
denote rectangular cartesian coordinates in a plane perpendicular to the 
axis of the wave-guide. Then & satisfies the equation 


: += + ye = %, (35) 


where y? k?— B2, 


Let x+iy = g(€+in) be some conformal transformation. Then 


25 arp 
O-w Cc se A > 
<2 4 a > | xyes — 0, (36) 
a C n° 
where a? = g’(€+in)g'(E—in) > 0. 


hen on the curves in the (£, 7) plane which correspond to the wave-guide 
walls % will obey the same conditions (i.e. 4 = 0 or &//én = 0 as appro- 
priate) as on the walls themselves. 
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Multiplying by % and integrating over the cross-section © in the (é, n) 
plane . P (ene ay 


- x7? di = — a <a pb dx 
: 7g yoo eer 


-f (C2) +] 2 foe 


. 


~~ 





where I is the curve in the (€, ») plane which corresponds to the guide walls 
in the (x, y) plane. 
Either ¢% or &/én vanishes and so 


e 


[ [(Gb/eé)?+ (Gb/en)?] dd 
a : . (37) 
xf? dd 





This may be compared with the usual formula for cut-off frequency (11). 

It may be shown in exactly the same way as previously that this is a 
variational expression for y*. It follows therefore that it is possible to find 
the propagation constant for an arbitrarily shaped guide provided that 
either % or é/én vanishes on the walls. The first condition is true when 
ys = E_ irrespective of any conditions on the nature of the field. The second 
condition applies when the field is a TE field and y; = H,. These conditions 
between them cover all possible kinds of field, as modes may be divided into 
four categories, namely (1) TEM,t+ (2) TE, (3) TM, (4) mixed. For a mode 
of th. ‘rst type 8? = k? (12), for a mode of the second type H, exists, whilst 
for both the third and fourth modes E, exists. Thus all possible types of 
mode are covered. 

If it is possible to find a conformal transformation such that the guide 
cross-section may be expressed in the form 0 < € < &, 0 < » < mp, com- 
plete sets of functions in the €, 7 plane suitable according as & or éy/én 
vanish on the boundary are 


sinn7€/Ey, sin s7n/ no. 

and cos n7€/Ey, COs smn/ no. 
The critical frequencies may also be calculated in the manner indicated 
in section 7. 
9. Discussion 

It will be observed that there are only a finite number N of variables 
available (the b,). Were there an infinite number it would be possible to find 
the field distribution exactly. As, however, it is not possible, except under 
very rare circumstances, to solve an infinite number of equations in an 


+ A TEM mode is a mode of propagation such that the longitudinal components of both 
electric field and magnetic field are zero. 
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infinite number of unknowns, it is necessary to use only a finite number V 
and the solution will be a compromise between small N, which involves 
comparatively simple numerical manipulation and comparatively low 
accuracy, and large NV, which is more accurate but involves more trouble- 
some calculations. Fortunately in practical cases no harm is done in erring 
on the side of simplicity because the electromagnetic constants of a material 
are rarely known to a high degree of accuracy. For example a nominal 
dielectric constant of 2-5 usually implies that, if a piece of material be picked 
up at random, the dielectric constant lies between 2-45 and 2-55. Under 
such circumstances an accuracy in calculation better than 1 per cent. is 
meaningless and a satisfactory answer may be had by the use of about three 
terms (see Appendix). The underlying reason for this is that properties 
such as propagation constant do not depend upon the details of the physical 
structure and so an approximation is good enough. 

It may be remarked that the method just indicated applies equally well 
for propagation of sound waves along a pipe which is filled with various 
media. The boundary condition to be applied in this case is that the particle 
velocity normal to the pipe vanishes. This effectively means that a/én 
vanishes over the walls, where y¢ is the velocity potential. 


APPENDIX 


The following example, extracted from (13), shows the degree of approximation 
involved. 

If a wave-guide 0 y a is filled with dielectric (€, 4) in 0 < y < cand dielectric 
€ > f4) In ¢ y a, the propagation constant may be calculated exactly in the form 
3/k, where k? = w*ue, and h 27 /Ay where Ag is the wavelength in the dielectric (€ , 1). 
Now E, must vanish at y = 0, y = a, and so a suitable way of writing the field is 

, 
E > A, sin(nzy/a) 
n=1 
and the appropriate partial approximations are 
N 
En > A, sin(nzy/a). 
n=1 
If now the following values for the constants involved are taken 
Ag = a, €/€g = 2-45, c= a/2, 
the exact value is given by 
B/k 1-3585, 
and the successive values of 8/k for the various approximations are 
N ] 2 3 


B/k 1-2145 1:3497 1-3546 


N l 2 3 


“%, error 10-6 0-7 0.3 


( 
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A FURTHER NOTE ON DUAL INTEGRAL EQUATIONS 
AND AN APPLICATION TO THE DIFFRACTION OF 
ELECTROMAGNETIC WAVES 


By C. J. TRANTER (Royal Military College of Science, Shrivenham) 
Received 18 March 1953] 


SUMMARY 
A previous solution of the dual integral equations 
x 


( Cau) t u)J,(pu) du : g(p) (O « p< 1) 


i f u)J,(pu) du 1) (p > 1) 

0 
where G(u), g(p) are given functions of the variables indicated and f(u) is to be found, 
isextended to cover cases in which the order v of the Bessel function kernel is not zero. 
As an example, the solution is applied to a problem in the diffraction of electro- 
magnetic waves by a plane slit and some criticisms are made of a recent discussion 


f this problem 

1. Introduction 

30UNDARY-VALUE problems in which different conditions hold over different 
parts of the same boundary can often be reduced to the solution of dual 


integral equations of the type 


( G(u)f(u)J,(pu) du =g(p) (0<p<—]l) (1) 


0 


x 


l fl uJ (pu) du 0 (p - 1) (2) 


0 
where G(w), g(p) are given functions of the variables indicated and f(u) is 
to be found. The special case in which vy = 0 is of importance in problems 
with axial symmetry and I have given a solution for this case in a previous 
paper (1). 

[have recently considered several physical problems which can be reduced 
to the solution of the dual equations (1) and (2) when the Bessel function 
kernel is of non-zero order. My previous analysis can be extended to cover 
such cases and this extension forms the subject of the first part of the present 
paper, 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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In a recent paper Groschwitz and Honl (2 a) have discussed the diffrac. 
tion of electromagnetic waves by a narrow plane slit. This problem can be 
reduced to the solution of a pair of integral equations of the type (1) and 
(2), and in the second part of this paper I have applied my solution to the 
diffraction problem. The results obtained show major differences from 
those given by Groschwitz and H6énl, and some criticisms are given of certain 
parts of their analysis. 

2. The solution of the dual integral equations 

The necessary analysis is a direct extension of that given previously (1), 
a brief account of which is given here. The reader should easily be able to 
fill in the details by reference to the previous paper (1). 

Taking re(v) > —1, ma positive integer (or zero) and k real and positive, 
the integral a 
I(v,m,k, p) = ( WK*T omsz(U)d,(pu) du (3) 

: 0 
converges and its value is given by Watson (3 a). For p > 1 the integral 
vanishes and f(u) automatically satisfies equation (2) if we take 


x 


f(a) = @-* > ay, Sy samrn(™), (4) 


m=0 
and the coefficients a,, have now to be chosen so that f(w) satisfies equation 
(1). 
In precisely the same way as in the previous paper we find the result, 
which is not without interest in itself, 


u “i, +2m rat) 
; 1 
T(v+m-+1) cea \k1g¢ 2 5 
a. 5 _ | pr+1(1— 92)k-1.F (k-Ly. y+]. o2)T(up)de, (3 
FALL Tomek | PPA Falk+v, w+, p*Wy(up)dp, (5) 
0 


where F,, (a, y,%) = oF\(—m, a+m; y; x) (6) 





is Jacobi’s polynomial (Magnus and Oberhettinger, 4). 
Substitution from (4) in (1), multiplication by 
p” +1(] — p?)* VF (k+y, v+1, p”), 
where v is a positive integer or zero, integration with respect to p between 
0 and 1, interchange of the order of integration, and use of (5) give 


x 


> A», [ G(u)w-* SF, omsr(U)I sonai(u) du = E(v,n,k), (7) 


m=0 . 6 
where, for brevity, we write 
1 
P(v-+n+1) , " ; 
rt Bie at OS v+1(] — p2)k-1F (k§-+v, y+], p?) dp. 
Qk-1] (y+ rere | g(p)p ( Pp ) A Vy p ) 


0 (3) 


E(v,n,k) = 
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Equation (7) with » a positive integer or zero gives a set of simultaneous 
equations for the determination of the coefficients a, in the series for f(u) 
given in (4). Since (see Watson, 3 b) 


( WS) oman (U)d, ,(u) du = 0 (m ~ n) 


0 : > ’ (9) 


= (2v+4n+ 2k)-? (m = 2) 


9 


these equations give expressions in closed form for a,, when G(u) = u?*-?, 


m 
If G(u) is such that k can be chosen so that G(w)—u*-? is fairly small, the 
simultaneous equations can be solved by an iterative process for the impor- 
tant special case g(p) = Ap”, where A is a constant, as follows. 


‘ 


Making use of (9), we can write equation (7) in the form 


a, > | ee (2v+4n+2k)E(v,n, k), (10) 

where ; 
L (2v+4n+2k) | fu?-**G(u)—Liu Wd, .onsn(W)Ipsonsn(u) du. (11) 
When g(p) Ap”, use of (8) and the orthogonality relations for the Jacobi 


polynomials (see Magnus and Oberhettinger, 4) gives 
E(v.n,k) = 0 whenn > 0, 
E(v, 0, k) 2-*T(v+1)A/T(k+v-+1), 


so that the right-hand side of equation (10) vanishes for n > 0 and is 


2-AD(v-+1)A/P(v+k) for n = 0. 


equations (10) is then 


The iterative solution of the set of 


yI-kKP( 1)A . 
yd d a . . , ) 
a, 7 [d,,- ( nt T ( n ( n T , (12) 

D(v+k) 
where 6 0 for n > 0, 6, = 1 and 
x x 
’ , - ’ we v ‘ 
( Lo, ? ( i > ae m? ( n > ( m? etc . ( 13) 
m=0 m=0 


3. An application to a diffraction problem 

In their paper on the diffraction of electromagnetic waves by a plane slit, 
Groschwitz and Honl (2 a) take the slit to lie in the (x, y)-plane and to be 
bounded by the lines x -b. They represent the incident electromagnetic 
wave, with the electric vector parallel to the y-axis, by uy = e’**, and they 
show that the electromagnetic field is given by 


x 


| F(a)e iK(axr+y2) doy 


—@ 


U(x, Zz) 
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in the upper half-plane z > 0, and by 
io 9) 
iia, z) — etkz__p—ixz + | F(a)et(or-y2) dx 
in the lower half-plane z < 0. Here a?+ +? = 1 and the weight function 
F(x) satisfies the dual integral equations 


1 on 
[ 1—a?)! F(«)cos(kéx) dx + 7 ( (a2—1)# F(a)cos(kéx) da = 4 (\é|< 5), 
i (14) 
( F(«)cos(k&x)dxa = 0 (\€| > Db). (15) 


0 
Groschwitz and Hoénl attempt a solution of these equations for the case 
of a narrow slit (6 small), but their method of solution seems to be open to 
criticism. The solution obtained when the analysis of section 2 is applied 
to equations (14) and (15), which are a special case of the more general 
equations (1) and (2), is given below. 


Writing u = Kba, p = &/b, e = kb, (16) 
f(w) = wF(u/e), (17) 
G(u) = —i(e®—u*)* when wu<e and G(u) = (u?—e?)! when u > «, 


(18) 
and noting that cos pu = (mpu/2)'J_,(pu), equations (14) and (15) become 


, . ie? " ; 
| G(u)f(u)I_,(pu) du = — Bap) (0<p< 1), (19) 
i f(u)J_y(pu) du = 0 (p> 1). (20) 


Here ¢ is small and G(u)—u?*-? is fairly small if we take k = 3/2. With this 
value of k and taking v = —}, equations (16), (17), and (4) give as the 


required solution 


, <= J5,.4(Kdx) 
= — : ee ——— , (21 } 
€ Am 


The coefficients a,, in (21) are given by equations (12) and (13). The 
latter require the calculation of the integral 
- 
Linn = (4n+2) | {u-2G(w)— Du Dana (W)Jonsa(u) du, 


0 
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which, by use of (18) and (9), is 
9 ;p S 99 
Mie: (4n T =)(- tP nt Qmn)—S mn: 22) 
where € 
) fr 95 
notion Fn | (e°—u*)iu-* Jy, 41 (W) Jo, 44(u) du, (23) 
2 2),,-2 9 
,) Q, | (u2@—e?)Pu-2 J, 4 1(U) Io, 44(u) du, (24) 
(14 ; . 
nd §,,,, is unity when m n and vanishes when m + n. 
15 In the calculation of P,, ,, and Q,,,,, powers of e above the second have been 
neglected. This gives the final solution for F(«) correct to e? and enables a 
e case comparison of the present result to be made with that given by Groschwitz 
en to and Honl 
plied From (23), writing 2 ex, 
-neral 1 
ld | (l—a?)ta-*J,,,, 4 (€x) Jp, (ex) dx. 
0 
(16 
" Use of the power series in ex for the product of the two Bessel functions 
si Watson, 3 c) shows that P,,,, is negligible except when m = n = 0, in 
€; which case 1 
18 of 9\} 9 indi 
Pio fe? | (1—2?)§ dz = 47’. (25) 
come . 
0 
0 The approximation for Q,,,, is more difficult. By using the integral 
representation (Watson, 3 d) for the product of two Bessel functions and 
then Bessel’s integral (Watson, 3 e) for the resulting Bessel function (which 
») . . 
“V) is of even order), equation (24) gives 
1 ' 7 har 
1 this a: ° 
» the Q —, | cos 2(m—n)6 dé | cos 2(m+n-+1)d dd » 
0 0 
” (u?—e?)! ' , - 
| ——, cos(2ucos@sind) du. (26) 
72 
; - : 
») € 
Denoting the integral in w in (26) by I, we have, writing w = ex and 
The » >A 3 
w 2e cos @sin d, 
* (~?—1) * COS wx “COS wx 
} I : cos wx dx > +. 2 ne dx. (27) 
; r J (a?—1)! J x*(x2*—1)! 
1 i 


Y 
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Using Watson (3/) 


x 


COS wa . 
[ Gi = —4rY,(w) = —{y+log 4o}+ 4w*{y— 1+ log 4w}—..., (28) 
where y here denotes Euler’s constant, 0-5772..., and we have neglected 


powers of w (and therefore «) above the second. Integrating twice with 
respect to w between 0 and w, we have 


x 


[ l—coswx 


6 





— 2. | 2 ant 5 
Pope TD t= —4w{y—$+log4w}-+.... (29) 
1 


The substitution « = sec 8 shows that 


4m 
r dx [ 
—_—__—_. = | cosS dB = 1, 
| x?(22?—1)* j B dp 
1 ) 
and substitution from (28), (29) in (27) then gives, correct to ¢?, 
I = —{y+1+log 4w}— }w*{y—2+log4w}+.... (30) 


Remembering that w = 2ecos@sin¢, the evaluation of the integral in ¢ 
in (26) depends on the definite integrals 


i7 477 
| cos 2(m+n-+1)¢ dd, cos 2(m+n-+ 1)¢ logsin d dd, 
0 0 

477 }7r 


| cos 2(m-+-n-+1)¢sin*¢ d¢, | cos 2(m-+-n-+ 1)¢ sin*¢ logsin ¢ dd. 

0 0 
These are either elementary or can be found from results given by Edwards 
(5a and b). Using these and (30) in (26), we find 





4a 
| 4 2 eos? 
in FT | 1— — cos 2(m—n)é dé, 
: (m+n-+1)7, 2(m+n)(m+n-+ 2) 
0 


when m + 0, n + 0, and 
$77 
Coe = ; | [1+ 4€? cos*6{y—§+ log(4e cos @)}] dé. (32) 
"9 
The final evaluation of Q,,,,, therefore depends on the definite integrals 
har har 
| cos 2(m—n)é dé, | cos 2(m—n)6 cos?6 dé, 
0 0 
47 ba 
| cos*é dé, | cos?@logeos 6 dé. 


0 0 
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The first three are again elementary and the last is easily found from a result 
given by Edwards (5a). With these results (31) and (32) give 





Qnn = 9 when m ~n and |m—n| ~1 > 
; I € : 
Grn = |__| whens > 0 
4n-+-2 l6n(n+ 1) = 
" (33) 
—e* | 
Grin = Cnn = 3 when n > 0 
32n(2n—1)(2n+-1) 
= Zsa ‘ 1 | 
Vo0 : 5 T sé LY $+ logte} J 


The use of (25) and (33) in (22) now gives, correct to order e?, 


Linn = 9 whenm ~n and |m—n| 41) 
e* 
a — ——___. when n > 0 
l6n(n-+-1) 
9 
e& 2 
Miigsitisn — whenn> 0 -. (34) 


16n(2n— 1) 


9 


2 
tint —____ when n > 0 
: l6n(2n+ 1) 
.. 2a. 3 ] 1,.-' 
Loo = de*{y— $+ log(te)— 47} J 
From (13) we now have 
7 3 : Fe 2 2 
C Le*{y— ?+ log(te)— iz}, C, = —te’, (35) 
and all other C’s are negligible to this order. From (12) withk = 3,v = —}, 
A ie?/(277)', we obtain finally 
Ay bie?[ 1 —1e*{y— 3+ log(fe)—427}], a, = — pie, (36) 


with all the other a’s negligible. 


4. Some remarks on the solution of Groschwitz and H6nl 
In the present notation the results given by Groschwitz and Honl (2 6) are 


“ S|} €* | 2 


ie4 10 
a, | l+ 2)? 
oe OT 


with the other a’s negligible. There are therefore major differences between 


_ 
— 
So 
_— 
«| 2 
| 


the two sets of results and these may well be explained as follows. 











TRANTER 


Using a rather different method of solution, Groschwitz and Hon require 
to evaluate the integral 


2 
| 1\3) , 
| _ 1 —_— “ Jon 1q(€ex)cos(Cax) da 
J | x") 
[ eth ama das Din satendoanibnd 
a a +... | J5,,.,(€x)cos(Cax) da, 
J [2a? 804 © sais 3 
1 
where € = exand |x| < 1. They retain only the first term in the expression 


in square brackets, whereas no such approximation has been made in 
section 3 above. This then will be one source of difference between the 
two sets of results. 

Groschwitz and Hénl then require to evaluate the integral 


rv [ly 


m 9 “2m+l1 
x 


(ex)cos(la) dx, 


ae 


and they do this by replacing the Bessel function by its integral repre- 
sentation 


for values of t between 0 and 7. 
Writing p = ecost+(, q = ecost—{, they find 
x ax = 2) x 
* Cos pa * Cos ga . f sin pa . f singa 
20 — I da + _ “. da +1 | . Ps dx +1 2 da 


x j x an - 
1 1 1 


me 


They then proceed to evaluate this for small p and q by integrating by 
parts and using the series for the cosine and sine integrals Ci(p), Ci(q), 
si(p), si(q) which then appear. They work, however, as if p and q remain 
always positive, and this is certainly not the case when ¢ ranges between 0) 
and z. That there must be an error in this part of their work is apparent 
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from equations (2.15), (2.15a) on page 314 of the paper by Groschwitz and 
Honl; these equations give the value of the real integral 


. fs . 
Ri} . Ji (ex)cos(Ca) da 
2 


1 


snot being wholly real. 
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CYLINDRICAL WAVE-GUIDES WITH IMPERFECTLY 
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[Received 30 April 1953] 


SUMMARY 

The effect, on the propagation of electromagnetic waves in cylindrical wave-guides 
of arbitrary cross-section, of imperfectly conducting walls is calculated by a perturba- 
tion method using approximate boundary conditions at the walls. The results show 
that, in cases where a particular type of degeneracy occurs in the ideal guide, the 
introduction of imperfectly conducting walls removes the degeneracy ; in such cases 
definite linear combinations of E and H modes are propagated with definite 
propagation constants. This type of degeneracy occurs in rectangular wave-guides. 
It follows that the power-loss method of calculating the attenuation constant, which 
assumes pure F or H mode propagation, is not correct in these cases. The perturba- 
tion method gives not only the attenuation constant for each combination of modes, 
but also the value of the phase velocity, and the corresponding field components. 


1. Introduction 

THE object of this paper is to investigate the propagation of electromagnetic 
waves in a cylindrical wave-guide with imperfectly conducting walls. The 
solutions for perfectly conducting walls are well known, but the only way 
in which the theory has been extended to include imperfectly conducting 
walls is in the calculation of the attenuation constant by the so-called power- 
loss method. By this method (see e.g. Stratton, 1), the power loss in the 
walls is calculated on the assumption that the field differs little from that 
in a perfectly conducting guide. The method gives no information about 
the field distribution, or about the change in the imaginary part of the 
propagation constant. In the present paper these effects are calculated by 
a perturbation method. It is found that the power-loss method leads to a 
valid result for the attenuation constant, except when the corresponding 
solution in the perfect guide is degenerate (i.e. such that two modes have 
the same phase velocity). When this is so, the effect of metallic walls is in 
certain cases to resolve the degeneracy (as would be anticipated from the 
general perturbation theory for eigenvalue problems). This results in the 
propagation of certain definite linear combinations of the degenerate modes, 
with a definite propagation constant associated with each combination. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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2. Propagation of electromagnetic waves in cylindrical wave- 
IN guides with imperfectly conducting walls 

For a wave propagated along the infinite z-axis of a cylindrical wave- 
L) guide of arbitrary cross-section, we consider fields which have a propagation 
factor exp(twt— yz) in all components, y being a complex propagation con- 
stant associated with the field distribution and with the shape of the guide. 
It is convenient to write Maxwell’s equations in terms of the transverse 
electric field E. The medium within the guide is taken to be free space, 
and rationalized m.k.s. units are used throughout. Maxwell’s equations 
then are (Marcuvitz, 2) 


iw 
yB = —— kAE+V(V.kKA E), 

; ce 
yuides 
turba- iwB, = V.KAE, 
a yE, = V.E, (1) 
cases where E., B, are the components in the z-direction of the electric field 
finite intensity and the magnetic flux density, E and B represent the components 


in the cross-section of these vectors, k is the unit vector in the z-direction, 





and c is the velocity of electromagnetic waves in free space. E must also 
10des, satisfy the equation 


nts. -V(V.E)+kAV(V.VA E) = pE, (2) 


where p is a constant defined by 


etl p = y?+w?/c?. (3) 
The Analogous equations for the field components within the metallic walls 
way | areobtained by substituting c*/u{«—(to/w)} for c? in the above expressions, 
aie where «. and « are the relative permeability and permittivity of the metal, 
wer and o is the conductivity. 

1 the The problem is to solve these equations, subject to conditions of con- 
that tinuity at the wall surfaces and the correct boundary conditions at infinity. 
bout For a circular cross-section this problem can be solved fairly easily; a few 
’ the other cases could perhaps be treated with difficulty. In the general case 
d by the problem cannot be solved explicitly; but for large conductivities it is 
toa | known that the effect of finite conductivity can be represented approxi- 
ding mately by imposing suitable boundary conditions on the field inside the 
lave guide at the metallic surfaces (provided that the metal is thick enough to 
is il give complete screening). These approximate boundary conditions are 
the Stratton, 1) 5.0 = of, 

pes } E, = —nB.s, 

" where » = 48w(1-++7), 5 is the skin depth 2!(wupyo)-? which is small for 
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large values of co, and s is a unit vector forming a right-handed system 
with the vectors n and k, in the order n, s,k. The unit vector n is that in 
the direction of the outward normal to the boundary S of the cross-section, 
In terms of E we may write these boundary conditions as 


E.s=”vV.kaE, 
tw 


v.E= nf SB. n+s.V(V.KAE)|. (4) 


The corresponding boundary conditions for perfectly conducting walls can 
be obtained by putting 7 = 0, and these conditions are 


E.s = 0, 
V.E = 0. (5) 


Equation (2) in conjunction with the relevant boundary conditions can be 
satisfied only by particular values p, of the expression {y?+(w?/c?)}. 


The solutions in the perfect guide are well known. There are two sets of 
solutions given (e.g. by Marcuvitz, 2) as 


E(t) = —V#4,, 

E? = —kA Vd, (r = 1, 2, 3,...), “) 
where (V2+ p\)d, = 0, % = Oon 8, 
and (V2+ p2))b. = 0, wg = O0on SV. 


E\” corresponds to an E mode of propagation, whilst E® corresponds to an 
H mode, the values p{!’ and p\”) are real, and in the two modes B is given 
as a multiple of k \ E by the expressions 





B = = \E (£E mode), 
ye 

B—=~-kAE (H mode). (7) 
iw 


The vector functions E), E\?) corresponding to different modes are ortho- 
gonal in the cross-section, and may be taken as normalized, so that we have 


| ( Ey). ED dS = Cus = J | E®).E (2) dS 
cross-section 


and | f EM .E® dS =0 (for p,q = 1, 2, 3....). 


These pr aiiitian may be derived from (2) and the boundary conditions (5). It 


may also be proved that the E® and E® form a complete set of vector 


functions in the cross-section. 
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The solution in the imperfectly conducting guide is to be found by a 
nerturbation method. The basis of such a method is the assumption that 
to each eigenvalue in a perfect wave-guide there corresponds an eigenvalue 
in the imperfectly conducting guide, the difference between corresponding 
eigenvalues being of order 7. 


A field E’ which is a solution of the problem for imperfectly conducting 


walls can be expressed as an orthogonal series in the E,. Thus 
E’ = > a4,E,, (8) 
r 
where a || E’. E, dS, and the summation is extended over both F and 


H modes. 
For any fields E, and E, satisfying (2) with arbitrary values p, and p, 
respectively of p, it is found immediately from (2) that 
op, || E,.E, dS Pa | | Ey. E, dS + 


A A 


({ V.SE, \k(V.kKA E,)—E, AkK(V.kA E,)— 
} 


_E,(V.E,)-+E,(V.E,)}d8. 


Hence, using the divergence theorem, 


(p,—p») | | Eq-E, aS 


q f(E,.s)(V.KA E,)—(E,.s)(V.k A E,)— 
_(E,.n)(V.E,)-+(E,-n)(V.E,)} ds. 


IfE,, and E, satisfy the boundary conditions (5), this gives the orthogonality 
relations for the modes in a perfectly conducting guide corresponding to 
distinct eigenvalues 

If only E, E. satisfies these conditions (5), and E, = E’ satisfies the 
conditions (4) for an imperfectly conducting guide, then we have 

OT oe P((V.KAE)(V.KAE 
p’—p,) | | E'.E,dS =n [P= )_ 
J | iw 

iw 


? RB nt s.VV.KAE’)|! ds. (9) 


Ww } 


(E,.n) 





( 


We now consider the mode in the imperfectly conducting guide for which 


has a value close to a particular one, pp, of the set of values p,. Then 


p 
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p'—p, is not small except when p, = py, and therefore, from equations 
(8) and (9), a 
a, = | | E’.E,dS = O(n) for p, + po. 
‘A 
It is possible that to the eigenvalue py there correspond several linearly 
independent fields E,, Eog,..., Eo,, (ie. the mode corresponding to p, is 
degenerate). Equation (8) may now be written as 


Om 


r m 
E’ = Z Vor E,,+0O(n), (10) 
r=1 
where ay, = {| E’. E,, dS. This result when applied to (9) gives (neglecting 
4 
terms of order 7?) 


, =< V.KA Ey, (V.KAE 
(P —Ppo)4or = 7 > ) don| = p AY BA Ber) _ 


A 


lw 


if 


r 


(Eyy-n)) Ey, n+ 2-s.V(V.k \ E,,) | ds, (11) 


where r’ = 1, 2...., m. 

Equation (11) gives m linear homogeneous simultaneous equations for the 
m constants d»,, in which p’ appears as a parameter. These equations have 
a solution other than a), = 0 (r = 1, 2,..., m) only for those values of p’ 
for which the determinant of the coefficients of the a), vanishes. If the m 
values of the p’ are distinct, there corresponds to each a set of the a), 
unique apart from an arbitrary constant factor; in such a case the modes 
are not degenerate in the imperfectly conducting guide. When, however, 
the equation for p’ has multiple roots, then the degeneracy is not com- 
pletely resolved in this approximation. To each n-fold root p’ correspond 
n linearly independent solutions, and the degeneracy of the system is 
n-fold. 

Now, to the first order in y, the propagation constant y associated with 
each value of p’ is given by (3), the corresponding fields are given by 
(10), (9), and (8), and all other field components can be derived from (1). 


3. The case in which the mode in the perfect guide is non-degenerate is 
especially simple. Equations (10) and (11) become 


KE’ = Ao, Eo,+ O(n) 
and 


V.KAE,,)(V.KA E, 
(p’—po) = nf , ui y 





lw 


Py 
8 


} V(V.KAE 
— (Ep -m)| (Ey) + $0 EA Ba 


9 I 
i tw 


| ds. 
J 
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Using (1) to express B and B, in terms of Ej), 


(p’— po) ” M) fiw B2+-y9(E,,.n)(B.s)} ds, 


where y) is the propagation constant in the perfect guide. For a frequency 


above cut off, w* > PoC"; so that y 


Yo = 18, with B, real. Equation (3) gives, 


to first order in », 


Ls P 
, re(y) 2B i1m(p — py). 





-0 
L . oe UF 
Hence Xx = 55 1M! » > {tw BE+-78,(E,,.n)(B.s)} ds}. 
-Vo . 
Using the results (7) 
\ a ) (B.s)?ds| for E modes, 
2B, Ww J 
and — im] inw > [B2+(B.s)?]ds| for H modes. 
=Po | : 


The definition of « by the power-loss method (see e.g. Stratton, 1) is 


4 re) f é/ #. n ds 


s 
\ . = 


{fe Hk dS) 
A 

where € and #& represent the complex 3-dimensional fields in a perfect 
guide, and # denotes the complex conjugate of A. Usingin the numerator 
the value of the tangential electric field obtained by the approximate 
boundary conditions (4), and applying also the orthogonality conditions 
and the results (7) to the transverse field components, we find for an H mode 


Inw [ ‘ | 
bre! B.)2-+(B.s)?] dsl, 
2 | Bo p K Ne 


and for an E mode 


tt 4,4 
5 Te\ 7 4) °”_(B.s)? ds}. 
' | a W } 
Thus, remembering that 7 }dw(1+7), 5 being a real quantity, the two 


methods lead to the same result in this case. 


4. The effect of degeneracy may be shown by considering the cases 
of circular and rectangular wave-guides. In a perfect circular guide, for 
all modes except those with rotational symmetry, there are two modes 
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corresponding to each eigenvalue (i.e. two polarizations). This effect arises 
from the circular symmetry of the problem, and this symmetry is not 
destroyed by the imperfect conductivity. It is therefore apparent that 
such a degeneracy will not be resolved in the circular guide with imper 
fectly conducting walls. The attenuation constant and the field components 
can be calculated from the formulae for a non-degenerate case. 

In rectangular guides, twofold degeneracy arises because E,,,, modes 


mn 


and H,,,,, modes have the same eigenvalue. We may write in this case, from 


(6), (10), (11), E’ = —a,Vd—a,k A Vb+ O(n), 








Ue ie r iw ap Od . Cus 
whence (p —po)a, = —7 , a, —+a,—) ds, 
c® On\ “en és 
Aa( ib) _ ob Me y ey 
(p’—po)az = nn — + y+ \ ds. (12) 
an iw cs ¢2 en iw Cs | 
In the rectangular sie of side a and b we have (Marcuvitz, 3) 
‘men? n2n?\-3 . max. nay 
Pnn = 2(ab) " xy + oy : sin — sin = 
a* b? a b a 
2,8 er (96, = BBs 
272 n*n?\-2 mre nm 
nn = 2(ab)-? — cos —— cos— Y 
a b? a b 


and the corresponding eigenvalue is 


9 9° 9 9 
mn  n*21" 
re ee a ee 
aa b? 


In the perfect guide above the cut-off frequency wy, when w > wy = Cp, 
1 ee 
Yo = U(w*—wy)*/c = to. 


The coefficients in (12) may be evaluated, the resulting equations being 


w?(n22a |) m2n?2h j w?mn7 q 
“i 
=(—,- + ——]—A]a,+ | ——— (a—b) Ja, = 0 
c2\ ob? a cab- 


and 
2mnr 202 272 
aaa Do [+ FE) rato se—alay = 0, (9) 
where hws tetbeus co(p'— Pe) 


4c? n 


A , 4c%1 —3 
Hence, from (3), y+ p= c ( 2 | pe d. 
we ab 24) wo 


It follows from (13) that, in the case of a square guide, pure E and H modes 
are propagated, with different propagation constants. 
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Equations (13) are satisfied by two values of the ratio a,/a, corresponding 
to two values of A. The expressions for these are involved and need not be 
given here. 

The propagation constant y x+i8 in the imperfect guide is given by 


3 Bj —re(p’ — po) + O(n") 


and \ a 
3 


im(p’ — py) + O(n?). 
0 
We are now able to calculate the field in the wave-guide and the associated 
propagation constant. 

The attenuation constants and the ratios a,/a, for the two modes are 
given as functions of frequency in the figures, for the case 


a 2b, m=—n= l. 


The attenuation constants for the ‘pure E’ and ‘pure H’ modes are shown 


lor comparison. 





| 





(a) fe} 
\ 
| 
(d) 
“aS \ 
a 
a 








5 678910 20 30 40 5060 80 100 wy, + 
J Wo 
Fic. | 
Curves (a) and (b) show the attenuation respectively of the combination 
modes given by curves (1) and (2), Fig. 2. Curves (c) and (d) show the 


respectivel. 


f H,, and £,, modes as calculated by the power- 
loss method. 
Another possibility of degeneracy arises in the rectangular guide for 
suitable ratios of side length. We may find for certain rational values 
\0 two different integer pairs (m,n) and (m’,n’) which give the same 
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propagation constant. In the square guide, for example, (m’,n’) = (m,n), 
It may be seen by the inspection of (13) that such pairs are not coupled, 
and therefore, as in the case of a circular guide, the two cases may be 
treated independently. 
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WHAT IS HAMILTON’S PRINCIPLE? 


By HAROLD JEFFREYS (St. John’s College, Cambridge) 


[Received 12 February 1953] 


SUMMARY 
In non-holonomic systems the so-called work function does not express the whole 
of the forces acting; the Hélder—Voss modification of Hamilton’s principle allows 
correctly for this fact. 


In a recent paper R. 8. Capon (1) quotes with approval an argument of 
Hertz purporting to show that Hamilton’s principle is not satisfied for non- 
holonomic systems and in particular does not yield Lagrange’s equations 
with undetermined multipliers. He also summarizes an argument of Hélder 
and Voss that does derive these equations from a different statement of the 
principle. This argument is given more fully by Whittaker (2)and by Jeffreys 
and Jeffreys (3). Capon condemns it on the ground that ‘no reasons, based 
on the theory of the calculus of variations, are given for the wide departures 
from the generally accepted formalism which are made’; the ‘generally 
accepted formalism’ being defined by the condition that ‘the comparison 
paths are required to satisfy the conditional equations, the displacements 
being free’. As full reasons for the departures are given in the works cited, 
it appears that Capon considers them unworthy of mention because they 
conc.n the physical conditions to be satisfied. I propose to show that 
Hertz’s argument fails to express these conditions and the Hélder—Voss 
argument succeeds. 

The fundamental point is that in the standard form of Hamilton’s 
principle the forces are supposed to possess a work function W, in the 
sense that however the system is moved from one position to another the 
forces do the same amount of work. In this strict sense no non-holonomic 
system has a work function, because in general the displacements would 
involve friction and the work done by the frictional forces would depend 
on the path. But if the paths are restricted by the condition that the 
frictional forces do no work, there may be a unique work function subject 
to this restriction. In a problem of rolling it is possible to consider paths 
where one body is lifted out of contact with the other and placed directly 
in its final position. It then becomes natural to ask what modification of 
Hamilton’s principle is required if the work function is defined in this 
restricted way. The need for some modification is plain. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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Lagrange’s variational principle leads tot 
hi 
| (87'+Q,89,) dt = 0 (1) 
to 
to the first order, for all variations of the positions of the system such that 
dq, = 0 at t = ty, t,. The generalized forces Q, refer to the actual motion. 
If part of the forces are derived from a work-function W, and the rest are 
frictional forces such that in a displacement 59, they do work R, 8q,, the 
principle becomes 
| (8(7-+W)+ R,8q,} dt = 0, (2) 


and since this is supposed to hold for all 5q, (subject to suitable conditions 
of differentiability) the usual methods lead to 


for all s. 

Suppose that the velocity of slip is a,(q)g,, which is to be 0 in the actual 
motion, and that the force in the direction of slip is R; then in a dis- 
placement involving slip the work done by R is Ra,5q,, whence 


R, = Ka, (4) 


and (3) are a set of equations with the undetermined multiplier R. 
We can also proceed as follows. Consider the displacements restricted 
so that . , 
: a,0q, = 0. (5) 
Then for any displacements satisfying these conditions the work done by 
all forces is the increase in W. But also in such displacements R does no 


(5 oL @al\, , 
d OL 8) 4 0 
dt Cd; Gs 


work, and we have 


subject to (5), and then by the method of undetermined multipliers we get 


d ol ol - 

— ine aN (7) 

dt @q, oq, 
Thus A is identical with R. This is the Hélder—Voss method. The reason 
for restricting the variations is that for other variations (5) would not hold 
and the work done by friction would have to appear in Q, 8q,. 


t+ The spmmation convention is used throughout. 
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In the Hertz method, on the other hand, the variational equation would 
be replaced by 
5 | (T+ W—nd)a,q,) dt = 0 (8) 


for some A(t) and all possible variations 6g, subject to dt = 0, dg, = 0 ar 
the termini. If in fact 


a, dq, dF (q), (9) 
we can write _— % 
6 | AF dt = —[A6F]+ | oF .X dt 
| Aa, 8g, dt, (10) 


which agrees with (2) for X R (instead of A = R as for the usual method 
of undetermined multipliers). Thus this method would work for a holonomic 
system, but in general for a non-holonomic system 


8 [ Aa,g, dt 4 [ Ra, 5q, dt, (11) 


and the principle in this form fails to express the dynamical conditions. 
Capon describes the Hélder—Voss argument as incorrect, but the only 
reason given against it is that it does not adopt the usual formalism. A 
sufficient reply should be that it attends to the physical principles and the 
usual formalism does not. The essential point is that the work function is 
defined in a way that represents only part of the forces; consequently there 
was never any reason to suppose that 6 { L dt = 0 for all variations of the 


path, even if the varied path satisfies the rolling condition. 
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VARIATION PRINCIPLES IN DYNAMICS 
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SUMMARY 
The original form of Hamilton’s principle refers specifically to holonomic systems, 
’ and various generalizations are possible when non-holonomic systems are considered, 
The paper discusses the validity of these generalizations, and their relations to other 
variation principles. 


1. IN a recent paper (1) R 8. Capon discusses the validity of Hamilton’s 
principle when the dynamical system considered is not holonomic. Capon 
asserts that Hélder’s theory is incorrect, and that Hamilton’s principle is 
not valid for non-holonomic systems. I do not think that Hélder’s theory 
is incorrect—indeed I regard Holder’s famous paper of 1896 (2) as one of 
the corner-stones in the theory of the variation principles—and I hope in 
the present paper to clarify the points at issue. 

The original enunciation of Hamilton’s principle refers specifically to 
holonomic systems. When we consider non-holonomic systems we can 
generalize Hamilton’s principle in (at least) two different ways: each of 
these generalizations is equivalent to the original principle when the system 
is holonomic. The first point to notice is that for non-holonomic systems 
one generalization gives a theorem that is true, and one gives a theorem 
that is false. 


2. Let us first consider holonomic systems. Hamilton’s principle can be 
expressed in two different forms. We consider an orbit of the system, i.¢. 
a dynamically possible path, from a configuration P, occupied at time ¢, to 
a configuration P, occupied at time ¢,. The termini and the instants of 
departure and arrival are fixed. We can picture the motion as the motion 
of a representative point in the n-dimensional configuration space, and s0 
long as we are dealing with a holonomic system we can take n equal to f, 
the number of degrees of freedom of the system. We can choose Lagrangian 
coordinates 4,, qo;---, Y,, and the orbit is a curve 


i 


Git) (r= 1, 2;..:5 n;t, <t< t,), 


where the functions q,(t) are of class C,,+ satisfy Lagrange’s equations of 


motion, and take the values corresponding to P, and P, at ¢, and f, respec: 
in a region R of the (2, Lgy-++5 Tm)” 


+ A function f(x, Xg,...,2,,) is said to be of class C,, 


space if it possesses continuous derivatives of order p in R 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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VARIATION PRINCIPLES IN DYNAMICS 339 
tively. We can enunciate Hamilton’s principle in either of two forms as 
follows 


(i) We can construct a varied path by considering at each instant a 


displacement 6q, (7 _ = n) from the point qg, (r = 1, 2,..., n) on the 
original orbit to the contemporaneous point q¢,+-4q, (r = 1, 2,..., n) on the 


varied path. The functions 5q,(t) ¢ C,. The varied path, like the original 
orbit, leaves P, at ¢, and arrives at P, at t,, so 6q, = 0 (r = 1, 2,..., n) at 
t, and f,. Hamilton’s principle can be expressed in the form 
te 
[ sLdt = 0. 
ty 
ii) We can consider the family K of all possible paths (i.e. geometrically 
possible paths) of class C, which leave P, at t, and reach P, at t,: the orbit 
ie. the dynamically possible path) belongs to the class K. We can express 
te 
Hamilton’s principle in the form that { LZ dt has a stationary value for 
the orbit in the class K. : 


For holonomic systems the two forms are equivalent. 

In each case the converse theorem is valid: the conditions are sufficient, 
not merely necessary, for an orbit. It is the converse theorems that we 
ippeal to when we use Hamilton’s principle as a practical tool to determine 
the motion. But in the present paper we shall be mainly concerned with 


the direct theorems, i.e. with necessary conditions for an orbit. 


3. We now consider a non-holonomic system. Let k denote the number 
of degrees of freedom, / the number of non-integrable constraints. These 
numbers & and / are invariants of the dynamical system, independent of 
the coordinates chosen to describe the configuration. 

The characteristic property of the non-holonomic system is that, whereas 
rom a given configuration a manifold of n dimensions, where n = k-+-1, is 
accessible by geometrically possible paths, a manifold of only k dimensions 
is accessible by dynamically possible paths. 

The least possible number of Lagrangian coordinates is n, and there are 


! equations of constraint 
> A,,dq,+A,dt=90 (r= 1, 2,..., l), (1) 


imposed on the possible displacements. The coefficients A,,, A, are functions 


s 


OF (9), Yo,---. g,3t) of class C,, and the system of Pfaffian equations (1) admits 


no integrable combination. The virtual displacements satisfy the equations 


> A,,89,=0 (r=1,2,..., 2). (2) 


1 
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Any point of the n-fold configuration space is accessible from any other by 
a geometrically possible path: but the manifold [ of configurations acces- 
sible from a given configuration in a dynamically possible motion has only 
k dimensions. 

It is therefore evident that at any rate the converse of the proposition (ii) 
of section 2 cannot be true for non-holonomic systems, since, if two con- 
figurations P, and P, are given, P, may not lie on the manifold T of con- 
figurations accessible from P, by dynamically possible paths. In this case 
the family AK of geometrically possible paths leaving P, at t, and reaching 
P, at t, exists, and in general contains a member satisfying the condition 

te 
8 | L dt = 0, but this member is not an orbit. 

t; 

This leaves open the question of what happens if P, is a configuration 

which is accessible from P, by a dynamically possible path. In this case 
te 

the orbit exists: does it satisfy in the class K the condition 6 ( Ldt = 0% 
ty 

We shall see that in fact it does not. The proposition (ii) is not valid for 

non-holonomic systems, even if P, lies on an orbit through P,. 


4. We can establish these properties by reference to the equations of 
motion and the equations of constraint, but it is simpler, and perhaps more 
striking, to deal with simple concrete examples. 


(i) A particle of unit mass moves in space in a given field of force, and 
is subject to the single condition 


zdx — dy 0. 3) 


The axes are fixed rectangular axes. The Pfaffian equation (3) is not 
integrable, and the. system is not holonomic. In this case k = 2, / = 1, 
n a7 

We first prove that any point of space is accessible from any other by 
a geometrically possible path. To fix the ideas take P, as the origin (more 
precisely, P, is the configuration of the system when the particle is at 0) 


and P, as (%5, Ya, 22). We consider the path 


y=fw), 2=/f'), (4) 


+ The example (i) is of wider significance than appears at first sight. Any system with 
3 


3 


s 


three coordinates and one non-integrable constraint w XY Ay, 49, 0, where the coetii- 
1 


cients A,, are functions of (q;, 92, q3) of class C,, is equivalent to the system (i). For, in 


virtue of Pfaff’s theorem (3), the Pfaffian form w can be reduced to the form z dx—dy by 
an appropriate change of variables. 
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er by where f(x) €¢ C3. The equation of constraint (3) is plainly satisfied, and we 
eces- have only to choose f(x) with the properties 

only f(0) = 0. f'(0) = 0: Slxe) = He. I’ (Ze) = 2. (5) 

vn (ii That infinitely many such functions exist is evident by a simple geometrical 

easel wgument: alternatively, we can choose as a concrete example of such a 

- function one of the family 

5 Case f(a) = (3Y¥,—2_ X_)(x/a_)®—(2Yg—2_ Xo)(x/X_)? + Ax*(x,—2)?. (6) 

chi ii) A more practical example is provided by a rigid sphere, of radius a, 

iti rolling on a rough horizontal plane. In this case k = 3,1 = 2,n = 5. To 

prove that any configuration of the system is accessible from any other, we 

proceed as follows. Suppose that in the initial configuration P, a point A 

— on the surface of the sphere is in contact with a point A’ of the plane, and 

nee that in the final configuration P, a point B of the sphere is in contact 

a with a point B’ of the plane. Consider a great-circle are AB, denote its 

ength by aa, and let the length of the line A’ B’ be a(2nz-+8), where vn is 

at in integer and 0 < 8 < 27. On the sphere let the are AB be produced 

if necessary) to C, where AC is of length a8, and let D be the midpoint of 

the are BC. To reach P, from P, we first rotate the sphere about the 

ns of diameter through A’ until AB lies in the vertical plane through A’ B’. 

more | Wenext roll the sphere towards B’ through an angle 2n7+ }(a+-8), so that 

| D comes into contact with a point D’ of the plane. We next rotate the 

nd | Sphere through an angle z about the diameter through D’, and then roll 

~ | through a further angle }(8—«) until B comes into contact with B’. Finally, 

| we rotate the sphere about the diameter through Bb’ until we reach the 

| desired configuration P,. The five operations can of course be carried out 

)  inafinite time, the coordinates during the motion having continuous second 

pa | lerivatives with respect to ¢. (A line can be turned through an angle y in 

the finite time-interval from t = 0 to t = 7 by taking its angular velocity 

_= it time ¢ to be 30yf?(7—1?)*/7°, and the angular velocity and acceleration 
‘ . vanish both at ¢ QO and att T.) 

at 0) | 5. The manifold of configurations accessible from a given configuration 


geometrically possible paths has n dimensions, but the manifold T 
ccessible by dynamically possible paths has only k dimensions. 
Consider as a concrete example the system (i) of section 4, and take the 


simplest problem, that of free motion: the particle is subject to no forces 


ther thar 


1e force of constraint. The equations of motion are 


i Az, y A, Z 0, 22 = ¥, (7) 





(5) 
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If we project from the origin O at t = 0 with velocity (wu, 0, w), the solution 
can be expressed conveniently (if w ~ 0) in the form 


cosh dé—1 ; 
r= 0- ¢ . y = 6- are Sal z = sinhd, (9) 
sinhd sinhd 
where 6 = wt, sinh¢d = ut. If w = 0, there is uniform motion along Or, 


The points accessible from O by dynamically possible paths lie on the ruled 
surface I’ given by (9), and we notice that any point of this surface can be 
reached at the given instant ¢, (> 0) by suitable choice of the velocity of 
projection. The configurations accessible from a given configuration by 
dynamically possible paths lie on a manifold of k dimensions. 

A similar result holds for the system (ii). We take fixed axes with 0 in 
the plane and Oz vertical, and we choose as the five coordinates defining 
the configuration the coordinates (&,7) of the point of contact with the 
plane, and the Eulerian angles (6, 4,4) defining the orientation, relative to 
the fixed axes, of a trihedral fixed in the sphere. If we again consider the 
simple problem where there are no forces except the force of constraint, 
the motion is determined by the simple rule that the angular velocity 
(w1,@,w3) Of the sphere, relative to the fixed axes, remains constant. 
The Lagrangian coordinates (£,7,0,¢,%) are uniform functions of the 
three variables (w,t,wat,wst), and the manifold [. of configurations 
accessible from a given configuration by dynamically possible paths has 
three dimensions. 


6. Consider now in the light of the preceding discussion the form (ii) of 
Hamilton’s principle. The principle cannot be valid if P, does not lie on 
the k-fold [ accessible from P, by dynamically possible paths. Suppose, 


te 
however, that P, does lie on T°: is it true that 6 | L dt 0 for the orbit in 


ty 
the class of geometrically possible paths? The answer is in the negative. 
To prove this it will suffice to discuss the system (i) of section 4, and to 
suppose again that there are no forces except the force of constraint. There 
exists an orbit leaving O at t = 0 and reaching P, at t = ¢,. The first 
question we have to ask is whether this orbit satisfies the condition 
te 
8 | Ldt = 0. Now a geometrically possible path satisfying the condition 
f; 


te 
5 | Ldt = 0, and also satisfying the prescribed terminal conditions, 
ty 
satisfies the multiplier rule (4), i.e. it satisfies the Lagrangian equations 
derived from the Lagrangian function 


4 (#2 +97? + 32) —p(ze—y), (10) 
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where the multiplier » is a function of t to be determined. The equations 
are — F = * ° ~— ° 
i pez pz, y fh, Zz — px, 2 = ¥, (11) 
and there is a unique solution with the conditions 
L,Y, 2) (0,0,0), (x, ¥, z) (u,0,w), pe A, (12) 
itt = 0, where A is arbitrary. The orbit (9) does not satisfy the equations 
il) for any choice of A. This is easily verified: alternatively we may 
notice that the relation between z and ¢ derived from (11) and (12) is 
22 fw2*— 2uAz—(A2—u?—w?)z*)/(1+-2?), (13) 
which is inconsistent with (9), whatever the value of A, except in the trivial 
casew = 0. The orbit leaving O at t = 0 and reaching P, at t = t, does not 
nder | L dt stationary in the class of geometrically possible paths leaving 


0 att = 0 and reaching P, at t = t,. Hamilton’s principle in the form (ii) 
of section 2 is not valid for non-holonomic systems. It was for this reason 
that Hertz (5) rejected Hamilton’s principle as a basis for analytical 


dynamics 


The form (i) of Hamilton’s principle is valid for holonomic and non- 
holonomic systems alike. 
We prove this from first principles: there is no need for any coordinate 
system more elaborate than the Cartesian coordinates of the individual 
particles of the system. The fundamental equation of dynamics (6) can 


be written in the form 


> (m,#,—X,) dx, (0). (14) 


Here x, Vos+. r 


of the system, the mass of the rth particle is denoted indifferently by 


, are the V (= 3v) Cartesian coordinates of the v particles 
M39, Mg,_1, Mg,, aNd (2,, d2o,..., 5a) is an arbitrary virtual displacement. 
The importance of this form, as was first noticed by Lagrange, is that only 
the given forces X,, Xo,..., Xy appear in the equation; the forces of con- 
straint do not appear. The condition is both necessary for the motion 

id sufficient to determine it. 

We consider an orbit (i.e. a dynamically possible path) and build up a 
new path by taking at each instant a virtual displacement (52,, 529,..., xy) 
trom the point (a, 25,...,%y) of the orbit. The variations are contempora- 
neous, i.e. the varied configuration x+6x in the new path is occupied at 
the same instant ¢t as the configuration x in the original orbit. The variation 


0X is a virtual displacement at time t, arbitrary save for the conditions 
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that 6x vanishes at ¢, and at ¢,, and that the components 62), 529,..., dz, 
are functions of ¢ of class C,. Since the variations are contemporaneous, 








: d . . 
6, = —ix, (¢ = I, 2...., N) (15) 
dt 
and 
te te 
Pe N - oN ale aN 
| (87+ > X,d2,) dt = D, (m, a, dt, +X, dx,) dt 
* r . r 
ty i 
oN, 1 
A a. : 
=i > (m i, — 64, +X,dx,) dt 
oe dt 
i, r=1 
te 
N te PN - 
= > m,iz,bz,| — | > (m,#, ~X,) 82,] dt 
r=1 t a r=1 . 
1 
4 (16) 
The integrated terms vanish, since 5x, = 0 (r = 1, 2,..., N) at t, and t, 


and the integrand in the last integral vanishes at each instant in virtue 
of (14). Thus 

. V 

| (67+ ¥ X,dx,) dt = 0. (17) 


ti 


If the given forces are functions of position only, and are derived from a 
N 
conservative field, the Pfaffian form > X,dx, is a perfect differentialt—aV, 


r=1 
and we can write the result (17) in the more familiar form 


ty 
| dL dt = 0, (18) 
i; 

where as usual L = 7’—V. The proof given does not involve the question 


whether the system is holonomic or not, and Hamilton’s principle in the 
form (i) is valid for holonomic and non-holonomic systems alike. 


8. It must be observed, however, that if the system is not holonomic 
the varied path described in section 7 is not in general a geometrically 
possible path. Again a simple example will serve to prove the point, an 
example which was given by Holder (2) and has probably been rediscovered 

+ Actually it is not necessary to assume that }) X, 62, is a perfect differential for arbitrary 
values of 52,, 52 ,..., d2y, but only for values of i 52xq,..., dy Corresponding to a virtual 


N 


displacement. It would suffice, for example, if }} X, 62, could be expressed in the form 
r=} 


—3V+Aw, where w 0 is a Pfaffian equation of constraint. 
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by every worker in this field since his time. We consider a particle moving 
in space and subject to the constraint 
w =adx+bdy+cdz= 0, (19) 
where a, 6, c are functions of (x,y,z) of class C,, and the Pfaffian form w 
joes not admit an integrating factor. The original orbit surely satisfies the 
condition (19), and so by hypothesis do the variations from it, so we have 
az+by+cz = 0, (20) 
10x -+ bday 1 ¢sz— 0, (21) 
ind the question is, Does the varied path also satisfy the condition? Let 
is assume that it does, i.e. we assume 
o(ax+ by- Cz) 0. (22) 
Now (21) is satisfied at each instant, so 
d - - 
—(adx + bdy + dz) 0. (23) 
ul ; 
Subtracting (22) from (23), and using (15), we find 


de 


da . db oa 7 - a ” 
| Ox ey oY y 6b) 4 | Oz : 8c} — 0). (24) 
\dt dt : 7 dt 
If we write (24) in extenso. we obtain 
a cc F : ob - =? 
1/02 >Oy) | }(2 Ox — © OZ)+ | ~ (voy —©~ Ox) 0 
i Cz 2 CX CX 
(25) 
But from (20) and (21) 
1 Oz = 8y : 20x Oz x jy y ox (26) 
a h e 
so (setting aside the trivial case in which 6x/% = dy/y = 62/2 and the 


lisplacement is along the orbit itself) we find, from (25) and (26), 


af: . cn o( “) \ f= | 0. (: 
CY C2 Cz Cx Cx cy 


But (27) is false, since (19) is not integrable, and the assumption that the 


wied path is geometrically possible leads to a contradiction. 
\s a second illustration, consider the rolling sphere mentioned in section 


t (ii). The equations of constraint, the so-called ‘rolling conditions’, are 


di acosédé— asin@sinddd’s = 0, (28) 
dy —asinddé + asin@cosddy 0. (29) 


Suppose that the original motion is one of rolling along Ox with no spin 


bout the vertical. If we consider a displacement in which at each instant 
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y = 0, the new path is possible! For the points of contact on the sphere lie 
on a circle of circumference 27a cos «, which is equal to 27a to the first order 
in a. But in general the varied motion is not possible without skidding, 
In the original motion 6 = ¢ = 37, » = 0, € = ay, and the variation js 


such that ebb & $n—a80 = 0. 


The equations of constraint (28) and (29) are satisfied in the varied path if 
5(€—a cos 46 — asin Osind yw) = 0, 
8(7—a sind 0 + asin# cosdy) = 0, 
which are equivalent (when we remember the values on the original orbit) to 
d 


“(86a 5h) = 0, — F(5n—a30)—apd$ = 0, 
a ( 


and these are not satisfied, except in the special case already mentioned 
when 6¢ = 0. 

Thus the varied path contemplated in section 7 is not in general a geo- 
metrically possible path. There was no reason to suppose that it would be 
possible. But this fact in no way invalidates the principle established in 
section 7. 


9. So far we have been concerned entirely with contemporaneous varia- 
tions. We now consider the situation that arises when a point P occupied 
at time ¢ on the original orbit is correlated with a point P’ occupied at time 
t+-dt on the varied path. 

Another generalization is also desirable at this point. The examples 
discussed hitherto have all been catastatic systems, i.e. systems in which 
the coefficients A, in the equations of constraint (1) are identically zero. 
We now abandon this restriction and consider the general case of acatastatic 
systems, in which the coefficients A, are not identically zero. As a trivial 
example of an acatastatic system we can consider a particle moving on the 
smooth floor of a lift which is rising with prescribed speed W (= W(t). 
Taking fixed rectangular axes, with Oz vertically upwards, the equation 


of constraint is dz— Wdt = 0. (30) 


This system is holonomic. As an example of a non-holonomic acatastatic 
system, suppose that the rough table on which the sphere rolls in example (ii) 
of section 4 is made to rotate about Oz with an inexorable angular velocity 
Q (= Q(t)). The equations of constraint, using the same axes and co- 
ordinates as before, are now 
dé —acos¢dé — asin @sindds + yn Qdt = 0, (31) 
dyn —asin¢ddé + asin@cosddys — €Qdt = 0. (32) 


These equations admit no integrable combination. 
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The characteristic properties of a catastatic system are (i) that ¢, = 0 
r= 1, 2,..., n) is a possible velocity, and (ii) that the class of virtual 
displacements (equation (2)) is identical with the class of possible displace- 
ments (equation (1)). Foran acatastatic system the two classes are different, 
and indeed, in general, they have no member in common. In the trivial 
example of the particle on the floor of a lift the possible displacements 
satisfy (30) and the virtual displacements satisfy 6z = 0. In other words, 
the possible velocities all have vertical component W, the virtual velocities 
all have vertical component zero. 

When we enunciate a variation principle which is to be valid for acata- 
static systems we must take care to specify in precise terms what kind of 
variation is involved. In Hélder’s principle (section 10) the variations are 
virtual displacements, although the variations are not contemporaneous. 
Thus the variation 5q is a member of the class of virtual displace- 
ments at time ¢, but the correlated peint q+6q on the varied path is 
occupied at the instant ¢+-d6t. The somewhat artificial nature of this 
set-up suggested to Voss (7) a variation principle in which the variation 5q 
isa possible displacement of the system in the interval 6¢: the correlated 
point q-+-8q on the varied path is occupied at time ¢+-d¢ (cf. section 11). 
[It seems that Voss himself failed to realize the importance of his discovery, 
and he left the problem with little more than a bare enunciation of a valid 
principle of this type. In fact a generalization of Hélder’s and Voss’s 
principles leads to one of the most fruitful and comprehensive of all the 


variation principles (section 12). 


10. To discuss the general case of an acatastatic system we introduce 
Lagrangian coordinates q,, Y»,---, ¢,, Where (= k+l) is the least possible 
number of Lagrangian coordinates. We have 

&, = fle d.nGs & = 1,4...., BD (33) 
where f, € C,, and our first task is the expression of the fundamental equation 
(14) in terms of the Lagrangian coordinates: we have to translate (14) from 
the language of x to the language of g. The operation is familiar, and leads 
to the equat ion 


yf = Fe | (34) 
a \dt\éq,, oq,} 


where (5q,, dgo,..., 5g,,) is an arbitrary virtual displacement: the components 
89, 3qp,.--, oq, Satisfy the equations (2). 

Consider then an orbit of the system, and the varied path in which 
the point q occupied at time ¢ in the original orbit is correlated with the 
point q+8q occupied at time ¢+8¢ in the varied path. The functions 
091, go,..., Og, ot € Cy. 
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We have 


ty ts 
- r Won, OL | 
[oat = | (Caer > Mag + Say.) at, (35 
© 7 \ ot — Cd, cq, 
ty ts 
where > denotes summation from 1 to n. Now 
ie @« — 
0g, = = 09,—9,— ot, 36 
> di: tr aT (36) 


so, substituting for 59, in (35) and integrating, we have 


to 
- =< OF .. |‘ (ol . oL | 
Laii= > tel a a — (> 42) oa t— 
a a ve Pe a Fad 


5L) 
- | Z aileq) — = 0q, dt. 


[f, therefore, the varied path has the same terminal points as the original 
orbit, so that 6q vanishes at ¢t, and t,, we have 


amd dy | 
it 
ty ; ty 


te 
[ sL+(5 q, 


ty 


ts 
= CLs \ ay _ ” os _ eL) ae 
N45 aot — ay Bt at - | Dk <3 ag Ode 


The integrand in the second member vanishes at each instant if (equation 
g | 
(34)) the variation (59,, 5qp,..., dg,,) is a virtual displacement at time ¢. In 
this case therefore 
to 


| sL+(> he 3 © 2 yl a = 0, (37) 
‘| — “ cq,) dt ct | 


. 


ty 


and this is Hélder’s principle.+ The crucial point is that the displacement 
5q is a virtual displacement from the original orbit at time ¢, but the point 
q-++4q is occupied in the varied path at time t+8f. We notice the contrast 
with Hamilton’s principle, where the two motions from P, to P, occupy 
the same interval of time. 

If the system is catastatic and if L does not contain ¢ explicitly, con- 
ditions which are fulfilled in most cases of practical interest, there is an 
integral E = h of the equations of motion, where E = ¥ q, (éL/éq,)—L. 
This is true whether the system is holonomic or not. Let us now consider 
the class of variations in which the condition E = h is satisfied on all paths 
of the family. £ is constant, not only on the original orbit, but also on the 


t The form (37) is somewhat more general than that originally given by Hélder. 
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varied paths, and indeed it has the same constant value as it had on the 


circumstances 


5L = (> a) 


riginal orbit. In these 


leads to 


or 
and |v! 


[a(S Ss. J+ (>. in) FT st! dt 0, 


i 


which is equivalent to 


sf (> 4 i, <) at 0. (38) 
cq, 
This is the principle of Least Action. 
11. We now return to acatastatic systems. If (8q,,5q9,...,89,) i8 a 
possible displacement in time df, and (q,,q.,...,g,,) is a possible velocity, 
i (89, —4, 3t, 8¢,—G2 t.,..., 84, —G,, 5t) 
sa virtual displacement. The proof is immediate, for 
SA 4,.@&@=@ (¢ = I, %...., §), (1) 
=1 
S A..g,.+A, 0 (r i an 1), (39) 
l 
whence > A,,(5¢,—G,5t) = 9 (r = 1, 2,..., l), (40) 


which proves the theorem in virtue of (2). The important case, and the one 


that we need for our immediate purpose, is that in which (4), Go,.-., q,,) is the 


velocity in the actual motion. 


Consider now an orbit of the system, and a varied path, the point q 


occupied at time ¢ in the original orbit being correlated with the point 


q+8q occupied at time t+-d¢t in the varied path. We calculate the corre- 


? have 


sponding 


change in | Ldt. We 
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and by (36) this becomes 


"ie OL. . oLd. oL d ..) 
iad} oo 
| a’ +2 04, Ir oss éq, dt 4r os Gn 04, —f dt ad dl 


te 
ij d(y els iS. OL as DEE ol), | 
| al> oq, ™_ (> fr og -1) a LX \dt 5a) | 4 


; . OL , OL) . 
+{5(> dt a 1) + et | | 
Now consider the coefficient of 5¢ in the integrand: we have 
¢ _ Ols Ld 4 OL, 
aD. tag —4)+ = D (alga) ae 
and using this result we find 


te ty 

7 a! fs d (eL ol 

8 | L dt = (> ad Bt) tee [ > (5 ( )- o"'\(5q,—4, dt) | de. 
i, h t; - 


r 
1 


r? 


anata | 4, eq,) 


If now 6q is a possible displacement in time 8¢, the integrand in the last 
integral vanishes at each instant in virtue of (40) and (34), and in this case 


ty th 

If we choose the variations so that 5q,, 5q,,..., 5g, 5¢ all vanish at ¢, and at 

t,, then te 

5 [ L dt = 0, (42) 
ty 

and this is equivalent to Voss’s theorem. If the system is holonomic, the 

result is‘equivalent to Hamilton’s principle. 


12. In the theorem just proved we require not only fixed termini, but 
also fixed times of departure and arrival. The requirement of fixed 
terminal points is unimportant, but the requirement of fixed terminal 
times is a serious embarrassment. However, in the case in which the 
system is catastatic and L does not involve ¢ explicitly, we can find a 
generalization of the theorem in which the terminal times are not fixed 
but may be varied at will. The generalization can also be deduced from 
Hélder’s principle. In this case there is an integral E = h of the equations 
of motion, and if, therefore, we fix the terminal points, but not the terminal 
times, we have from (41) 

te ite 
5 | (L+h) dt = (h—E)8 | = 


t; 





te | ‘ ite | 
5 (tat —(> oH 54, — Bst\| 41) | 
5 L dt = (> dy 4 (41) 
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In this equation 4q is a possible displacement in time df, but here (since 
the system is catastatic) the possible displacements are the same as the 
virtual displacements. 
We thus arrive at the theorem mentioned 
te 
6 | (L+h) dt = 0. (43) 
4 

The conditions required for the truth of this result are as follows: the 
system is catastatic, L does not contain t explicitly, the path is a dynamic- 
ally possible path, 5q is a virtual displacement at time ¢t, though the 
variations are not contemporaneous, and, finally, the end-points (but not 
the times of departure and arrival) are fixed. I do not know when this 
theorem was first discovered. I found the result some twenty years ago, 
but I expect it has been discovered and rediscovered many times. It is a 
powerful result, and one of the most fruitful of all the variation principles. 
For example, the principle of Least Action follows immediately: if I does 
not contain ¢ explicitly, and we consider the family of paths on which E 
is constant and has the same constant value h/ as it has on the original 
orbit, (43) leads at once to (38). Another deduction is concerned with the 
change of the independent variable in the equations of motion from t¢ to @, 
where dt = ud@, u being a given function of q. But the detailed dis- 
cussion of the consequences of the theorem (43) lies outside the scope of 


the present paper. 
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PROJECTIVE RELATIONS IN PLANE KINEMATICS 


By W. J. DUNCAN 


(Department of Aeronautics and Fluid Mechanics, The University, Glasgow 
[Received 21 May 1953] 


SUMMARY 
It is shown that any descriptive property of the instantaneous centres and 
normals to point paths of a set of rigid bodies in coplanar motion is preserved when 
the system is projected on a plane. A similar result holds for motion on a spherical 


surface. Proofs of kinematic properties may be simplified by use of this principle. 


1. Introduction 

At first sight it would appear that projection can have no useful application 
to rigid plane kinematics, for the projection of a moving rigid body does not 
move asarigid body, except in trivial cases. However, as will be shown, there 
are useful applications of projection to the kinematics of a set of bodies in 
any single configuration. This depends on the fact that descriptive theorems 
concerning instantaneous centres and normals to paths remain valid when 
the plane kinematic system is projected on a plane. It may then be possible 
to project the given kinematic system into another of simple geometric type 
whose properties can easily be established. 

The proposition which is proved below can be stated formally as follows: 
Any descriptive theorem concerning the instantaneous centres and the 
normals to point paths of a plane kinematic system remains valid when the 
whole configuration is projected on a plane. 

By a descriptive theorem we mean, in conformity with usage in pure 
geometry, one solely concerned with the concurrence of lines and the col- 
linearity of points. It is to be understood that the kinematic system consists 
of rigid bodies which all lie and move in a common plane. An analogous 
proposition holds for motion on the surface of a sphere. 

It is shown here incidentally that important propositions of pure 
geometry, such as those of Desargues and of Menelaus, arise naturally 
in kinematic investigations and are readily proved by kinematic arguments. 
2. Preliminary discussion 

We consider » coplanar rigid planes in relative motion and identify each 
plane by one of the integers 0, 1,... (n—1). As a matter of convenience the 
plane 0 will be regarded as a reference plane. The instantaneous centre for 
the planes r and s is the point where the relative velocity of these planes 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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PROJECTIVE RELATIONS IN PLANE KINEMATICS 353 
is zero at the instant considered; this will be designated the point rs where 
ICS the order of the integers is immaterial (rs = sr). 
It is clear that we may assign arbitrarily the positions of the (n—1) 
instantaneous centres 01, 02,..., 0(n—1) and the (n—1) angular velocities 
| Jo» M290-++> P%n—a9 Of the planes 1...., (n—1) relative to the reference plane 0. 
—_— However, the unit of time can be chosen so that a selected angular velocity 
has the value unity and the number of independently assignable angular 
velocities is then only (n—2).+ The total number of arbitrarily assignable 
quantities on which the configuration of the instantaneous centres depends 
is N 2(n—1)+(n—2) 
1 wher 3n—4 (2.1) 
: since each centre has two coordinates. Whatever values may be assigned 
to these .V quantities, there can be no incompatibility of the motions of the 
planes. The angular velocity of the plane s relative to the plane r is 
We, Wag — Wo (2.2) 
thes and Wye o,,. (2.3) 
lies it Next let us consider the set of three planes 0, 7, s. The centre rs may, 
orems so far as its velocity is concerned, be regarded as belonging to either r or s. 
wher Considered as a point of 7, the normal to its path 
ssible on the reference plane is the join of rs and Or, 
c type while, considered as a point of s, the normal to 
its path is the join of rs and 0s. Hence the three 
llows points Or, rs, and s0 are collinear and, since any rt 
id the of the planes may be regarded as the reference Fic. 1. Instantaneous 
en the plane, we deduce the well-known proposition RE ae ee ee 
that rs, st, and tr are collinear. We shall describe this property as ‘collinearity 
1 pure insets of three’. The identity of the magnitudes of the velocity of rs relative 
eC to the plane ¢ when considered as belonging to r or to s leads to the relation 
sists see Fig, | 
7 w,_t) 78 We ts. sr, (2.4) 
- which is of a symmetrical form; here tr.rs and ts.sr represent directed 
Bi segments 
i Z Suppose now that we have four planes 0, 1, 2, 3. We may assign the 
' centres 01, 02, and then 12 can be assigned as any point on the join of 01 
ind 02 since an angular velocity ratio a 9:@,) corresponding to this can 
eat ve lound, by equation (2.4). Similarly, we may assign 03 and then take 
ce the 13 as any point on the join of 01 and 03 (see Fig. 2). The centre 23 is now 
re Io! lo put the matter in another way, the configuration of the instantaneous centres of 
lanes the wi 8 n dey is 2) ratios of angular velocities. 
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fixed as the meet of the joins of 02, 03 and of 12, 23. Thus all 6 centres of 
the 4 bodies are now determined. Moreover, if we take 6 points which are 
collinear in sets of 3, we see that they are a possible configuration of the 
centres of 4 bodies. By equation (2.4) we have 


Woo? W19 01.12:92.21 
AD: BD 

W39 + Wo 02.23:03.32 
BF:CF 


W190: Wao 03.31:01.13 


CE: AE. 


A 
01 


Fic. 2. Instantaneous centres of four bodies. 
Hence we get by multiplication 


AD.BF.CE _ , 
BD.CF.AE — 


AD.BF.CE 


- ——— > = coast ] ‘ (2.5) 
DB.FC.EA 


or 
This is the theorem of Menelaus, which has here been proved by a kinematic 
argument. 

Next take five planes 0, 1, 2, 3, 4. We can now assign the additional 
centres 04 and 14, the latter being collinear with 01 and 04 (see Fig. 3). 
The centres 24 and 34 can then be determined as intersections, as shown 
by the dashed lines. But 24 and 34 must be collinear with 23 (see the 
chain-dotted line in Fig. 3). This amounts to a proof of Desargues’s theorem, 
for the triangles 02, 03, 04 and 12, 13, 14 are in perspective, with 01 as 
centre. The points 23, 34, 24 are the meets of corresponding sides of the 
triangles and we have shown that they are collinear, in accordance with 
the theorem. 

In Fig. 3 the pairs of points (02, 12), (03, 13), and (04, 14) have arbitrary 


positions on arbitrary lines through the arbitrary point 01 and the remaining 
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points are determined as intersections. Hence the configuration is the most 
general which is consistent with the 10 points being collinear in sets of 3. 
We can generalize this and state the 
Proposition ; 

Any set of 4n(n—1) coplanar points which are ‘collinear in sets of 
three’ is a possible configuration of the instantaneous centres of n 


coplanal rigid bodies. 





Fic. 3. Instantaneous centres of five bodies. 


The propositions of this section all remain valid for motion on the surface 
fasphere, but the sines of arcs must be used in place of the arcs themselves 
n the equations for the angular velocity ratios and in the statement of the 


heorem of Menelaus. Also ‘great circle’ replaces ‘straight line’. 


3. The general projective property 


1 


The normal to the path on the plane r of a point P of the plane s is the 
oin of P to the instantaneous centre rs. Suppose then that we have a plane 
kinematic system of n planes, whose $n(m—1) instantaneous centres are 
necessarily collinear in sets of three, also certain points P, Q, R,... and the 
normals to their paths on certain planes. Let the whole system, in one of 
ts configurations, be projected on a plane. The projection of the set of 


instantaneous centres consists of }n(m—1) points which are collinear in sets 


e possible centres for n planes, by the proposition 


f three nence they a 
established in section 2. Also the projection of the normal to the path of 


ny point P is projected into the corresponding normal. Now any ‘descrip- 


tive’ property of a figure is preserved in projection. Hence any descriptive 
property of the instantaneous centres and path normals of a plane kinematic 
system is preserved in its projection. In view of the remarks at the end of 
section 2, this proposition is also valid for motion on the surface of a sphere. 
As an « xample we may take the following property of an articulated 


hombus ABCD, which is easily proved by elementary methods: if the 
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normals to the paths of A and C meet on BD, then the normals to the paths 
of B and D meet on AC. This property is descriptive and it therefore 
remains true for a general articulated quadrilateral (1), for such a quadri- 


lateral can always be projected into a rhombus. 


REFERENCE 
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SEXTIC EQUATION 

By G. MILLINGTON (Marconi’s Wireless Telegraph Company Limited) 
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SUMMARY 


\n iterative method is given for solving the sextic equation with real coefficients. 

By accepting in the final stage the task of solving by a direct method a cubic 

ation that may have complex coefficients, it is possible to make the process 

an initial choice of only one quantity. This advantage, together with 

e achievement of a rapid convergence, enables the method to compare favourably 
er iterative solutions. 


1. Introduction 

Tue technical advances of the last decade have greatly increased the demand 
for the rapid and accurate location of the zeros of polynomials of high degree, 
is is shown by the number of papers appearing on this subject. An out- 
standing contribution has recently been made by Olver (1) in a detailed 
liscussion of general methods applicable to equations of any degree. Among 
the direct methods he considers the root-squaring method at length, but he 
onfines his discussion of indirect methods to cases where an approximate 

ilue of a zero has already been found by a direct method. 

He points out that for equations of degree less than six, say, it is often 
better to use a method especially applicable to such equations than to use 
general method. There is a wide field of technical problems where poly- 
omials of degree six or not much greater are concerned, and for which 

irious iterative methods have been proposed. A good account of some of 
these is given in a Government Report on Servo-Mechanisms (2), where 
stress is laid on obtaining simple equations for performing each stage of the 
iteration. It appears, however, that for polynomials of degree higher than 
lour it is necessary to choose initial trial values of two of the quantities upon 
which the solution depends, both of which have to be improved by iteration 
belore the zeros can be found to a desired accuracy. 

1 
sextic equation which is an iterative process requiring an initial trial value 


if 


e present paper deals with an indirect method for the solution of the 


only one quantity. This simplification is achieved at the expense of 
elaborating the numerical work at each stage of the iteration, and of 
ccepting finally the task of solving a cubic equation that may have com- 
plex coefhcients. At first sight this may be thought unpractical, but the 
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method can be made rapidly convergent and is found to compare favourably 
with others in speed and ease of handling. 

Some aspects of the method may be capable of generalization, and in 
particular it may be feasible to extend it to the equation with complex 
coefficients, but the discussion will be limited specifically to the sextic 
equation with real coefficients. It includes, however, the quintic equation 
regarded as a sextic equation that is known to have a zero root and at 
least one other real root. 


2. Description of the method 
The sextic equation is first reduced to the standard form 


9 


284 qzitbz3+-cz?+-dz+e 0 (1 


as this greatly simplifies the subsequent algebra and numerical work. 
Instead of the usual procedure of isolating a quadratic factor, the poly- 
nomial in (1) is resolved into a pair of cubic factors to give the equation 


(23 +- az? B,2-4 ¥1)(2* x2 +-B, 2+.) 0. (2) 


[t will be noted that, if (1) has at least two real roots, the resolution in (2) 
can always be performed to make all the new coefficients real, while, if all the 
roots are complex, it is possible in (2) to make the cubic factors conjugate 
complexes, and in particular a will then be wholly imaginary. This suggests 
that it is useful to replace y, and y, by two quantities / and m such that 


l+-m ; 
V1 (3) 


and y. . (4) 


Not only does this simplify the actual comparison of (1) and (2), but, in 
the case in which the cubic factors are conjugate complexes, / is wholly 
real and m is wholly imaginary. 

From (1) to (4) it now follows that 


B,+B.—o? a, (5) 
I—a(B,—B.) = b, (6) 
B,B.—am = c, (7) 
1(8,+8,\—m(B,—B,) = 2d, (3) 


l2?—m? - 
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Equations (5) and (6) give 


: at+oa? I—b 
en a” (7) 
ita? I—f 
and es en (11) 


whence, from (7), (8), and (9), by rearrangement, we find 


Se 


2d i. a e 
Pala m(I x ] 0, (13) 
[2 4e+-m?. (14) 
When the roots are all complex, it is convenient to write 
m am’, (15) 
x ix’, (16) 


so that m’ and a’ are real, and (10) to (14) become 


8, = (5) (17) 
iat (5) (18) 

Qn 
a oo 
pane m8) ‘» 
[2 = 4e—m”2. (21) 


[t will be seen from (17) and (18) that 8, and 8, are conjugate complexes 
as required, while in (19), (20), and (21) all the quantities are real and the 
equations are analogous in form to (12), (13), and (14). It is therefore 
sufficient to confine the description of the process to one of these sets. 

Considering the case in which m and «a are real, it will be seen that 
changing the sign of m is equivalent to changing the sign of « and leaving 
‘unaltered, in fact to interchanging the factors in (2). It is therefore possible 
to adopt the convention that m is positive. If e is negative, when the 
equation necessarily has at least two real roots, it follows from (14) that m 
must be > 2\e| in order that / may be real. 

The iterative process is then as follows: 

i) Choose a trial value of m:; (ii) find 7 from (14); (iii) use (13) to find a: 
and (iv) obtain a new value of m from (12). 

















360 G. MILLINGTON 





































This process has apparently two serious drawbacks: 
(i) The combination of (14) and (13) leads to six possible new values of m, 
This is, however, found in practice to be an asset, as in the prelimi- 
nary exploratory stage it enables the possibilities of a trial value of m 
to be quickly appreciated and directs attention to the most promising 

line of advance. 
(ii 


— 


Equation (13) implies the solution of a cubic equation at each stage 
of the iterative process. It will be noted that it is of the standard 


reduced type x84 3nx-+29 = 0. “a 


It is the author’s experience that the analytical solution of this equation 
can be put in a form that is convenient for rapid and accurate computation, 
In the preliminary stage the equation does not have to be solved to high 
accuracy, while an approximate root can be quickly improved by the 
iterative equation derived by Newton’s method, namely, 

tha 2 ey 
4 3(vz+p) 
or Ly = ¥(x3—q)/(aZ+p). (23) 

Once the main iterative process has been started and attention has been 
confined to one root of the cubic equation, it will be found in general that at 
each stage the required new value of « can be found directly from (23) by 
using the new values of p and q and a probable value of « estimated from 
the values obtained in the previous stages. Thus in practice in the iterative 
process the direct solution of (13) needs to be done only in the exploratory 
stage and to a modest order of accuracy. 

Assuming that the process has been made to converge satisfactorily, the 
remaining problem is to find the roots of the cubic factors in (2). When the 
coefficients are all real, they can be reduced to standard form and solved 
directly or by iteration if preferred. When the cubic factors are conjugate 
complexes, only one of them needs to be solved to give one each of the three 
conjugate complex pairs of roots of the sextic equation. 

The cubie equation with complex coefficients can be solved by a direct 
method, but as this fact may not be widely appreciated, an outline of the 
method is given here. It is first reduced to the standard form in (22), where 
p and q are now complex. The usual procedure can then be adopted by the 


substitution ‘ 
imte=. (24) 

t 
leading to 8 = —q-+r', (25) 


where r= q?+p'. (26) 
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Qne sign or the other can be used in (25), but if \q?| > |p*|, so that, from 
26), |r!| ~ |q|, it is preferable to avoid small differences by choosing the 
sign that makes tf? ~ —2q. The three roots of ¢® are then extracted 
ind when they are substituted in (24) they yield the three complex values 
ofr. The achievement of a rapid and accurate solution depends on adopting 
, systematic method of extracting the square and cube roots, and with a 
ittle experience the method is found to be quite practicable. It may also 
be noted that, in general, (23) can be used when an approximate value of 


, complex root has been found. 


3. The convergence of the iterative process 

It is convenient to replace m in (12) by n so that n calculated from 
14), (13), and (12) may be regarded as a function of m. Confining attention 
toa chosen branch of the n(m) curve, the true value of m required, say mo, 
will be given by the intersection of the curve with the line n = m. A value 
fn that differs greatly from the m value from which it has been derived 
loes not necessarily imply that the m value is far from mp, as it may happen 
that the n(m) curve is very steep in the neighbourhood of m = mp. 

In the interests of rapid convergence it is not desirable to use the n value 
isanew m value for the next stage of iteration, quite apart from the question 
f whether this would lead to a convergent process. Instead, having chosen 
nthe preliminary work two neighbouring values m, and m, of m which 
refer to a branch of the n(m) curve that appears likely to intersect n = m 
ta value of my near by, a new value mz, is defined by the intersection of 
the chord through the points (m,,,) and (m,,”,) with n = m. 

The value of m, may be computed from either of the equivalent forms 


Mo.— Ne - 
Ms m,-+-— : (27) 
7 a— ] 
No— Mm. 
D Mes No : . (28) 
: ™ a’ —] 
where F _ —— . (29) 
s’ Ms—™M, 


but it is desirable to use (27) when s > 1, and (28) when s <1, so that m, 
an be computed as a small correction to the best available approximation 
10m, it being assumed that (mp, v») lies nearer than (m,,n,) ton = m. 

The value m, can be said to be obtained from m, and m, by asimple chord, 
t by interpolation or extrapolation according as (m,,”,) and (m.,, 7.) are 
opposite sides or on the same side of n = m. If mz has been found by 
iterpolation, the convergence of the iterative process is greatly expedited 
ithe next value m, is found as the intersection of n = m with a chord 
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through (m3, 7”) parallel to the chord joining (m,,n,) and (m,,,). In this M 
case m, may be said to be obtained by a parallel chord. 

The value of m, is then given by il 

mm, = m+ 33 (30) 
s—l 
/ Ng—Ms if 
or mM, = "st , (31) “ 
where s and s’ are given by (29). The computation is thus simpler than foll 
performing another interpolation, as the slope s will have been computed tria 
already for the preceding interpolation. crit 

The next value m, is then obtained from m, and m, by using a simple i 
chord. If this is an interpolation, it is followed by a parallel chord con- val 
struction, but if it is an extrapolation it is first followed by an interpolation, imé 
and so on. This sequence is found to converge much more rapidly than wh 
the successive application of the simple chord construction, and it has the val 
great merit of approaching the limit m, from both sides. It is designed to nal 
employ only the simple linear forms in (27), (28), (30), and (31), since a an 
second-order approximation such as Newton’s method, though possible, ” 
would be very tedious where n(m) is found from equations involving a 
quadratic and cubic forms. | 

In principle such methods may be applied to an iterative process in which 
trial values of two quantities have to be chosen initially, but the added 
complexity tends to offset the gain in the rapidity of convergence. It is lie 
the essence of the present treatment that these methods for assuring a pr 
rapid convergence become practicable when dealing with the single quantity of 
Mg, and justify the use of the less simple type of iterative process contained fu 
in (12), (13), and (14). 

It is implied that the curvature of the n(m) curve in the range under 
consideration is not so great that an extrapolation using the simple chord 
construction may produce a degraded result. In this event the initial values 
m, and m, would have to be improved by a closer preliminary inspection 
before the routine process for converging on to m, is carried through. Such 
an initial difficulty is, however, common to most iterative processes, and 
is not a weakness of this method as such. 

4. A numerical example 
As an example the following equation will be solved: 
p®+21p>+ 210p4+ 1260p? + 4725p?+ 10395p+- 10395 = 0. 

This equation arises in connexion with the design of a maximally flat " 

delay network and its roots have been given independently by Thomson (3) : 


and Darlington (4). 
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Writing p = 2—3'5, 
the standard form of the equation is 
26 +. 26-2524+ 3523+ 177-18752? + 215-25z+-193-046875 = 0. 


It was known on physical grounds that the roots must all be complex, 
so that (19), (20), and (21) were used with m’ in (19) replaced by n’. It 
follows from (21) that m 2e*, i.e. < 27-8. It is inadvisable to choose a 
trial value of m’ that is close to this limit, as the calculations would become 
critically dependent on the small value of / that would result. 

An initial value of 20 was chosen for m’, giving 1 = +19-3. The positive 
value leads to a cubic equation having two complex roots with large 
imaginary parts and a real root a’ = 3-05. The latter leads to n’ = 32-35, 
which is not very promising, especially as it is greater than the limiting 
value of m’ above. The negative value for / gives three real values for a’, 
namely 1-19, 6-29, and —7-49, which yield for n’ the values —414, 18-1, 
and 7:83. Of these, the second seems the best to pomeee, and taking m’ = 18 
and confining attention to this branch of the n’(m’) curve, the values 
l 21-17, a’ }-237, and n’ 18-75 are obtained. 


It will be seen that the points 
(m, 20,4 18-1) and (mj 18,5 18-75) 


lie on opposite sides of n’ m’ and form an excellent start for an iterative 
process to find a value of mj that must lie not far from 18-5. The value 
of a’ is already sufficiently well known for (23) to be used to find all the 
further approximations. 


The results of the iteration can now be given in Table 1: 


TABLE | 


x’ n’ | om’ obtained i by 
7 | 

6°29 18-1 Pecinainers work 

6°237 18°75 Preliminary work 
6:2509 18°5965 Interpolation 
62514 18-5911] Parallel chord 
1859071 5364 6°251436 18590701 | Extrapolation 
650 6°251435 18°590703 Interpolation 


Care has been taken to include enough figures at each stage to ensure the 
inprovement sought. In this example the slope of the n’(m’) curve is of the 
order of 0-3, so that the forms in (28) and (31) were used. It will be seen 
that the convergence is remarkably rapid, four stages only beyond the 
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preliminary work sufficing to give m’, 1, and a’ correct to six places of 


decimals. It follows from (3), (15), (16), and (17) that 


Y¥, = —10-326825-+7 9-295351 
and B, = —6-415220+ 7 4-451270. 


When the corresponding cubic factor in (2) is resolved and the roots are 

reduced by 3-5, the roots of the original sextic equation are found to be 
—3-73571-Li 2-62627, 
— 4-24836-+7 0-86751, 
—2-51593-4i 4-49267. 

These values conform to those given by Thomson and by Darlington and 
take them to a slightly higher degree of accuracy. 

As a check on the accuracy of the solution, it is interesting to solve the 
same equation using another branch of the n’(m’) curve. From the known 
roots the other possible values of mo can be calculated, and very roughly 
they are 14-6, 23-2, and 25-5. Inspection shows that mj, ~ 14-6 lies ona 


very steep branch of the n’(m’) curve and that m, = 14-8 and mj = 145 
give n,; = 9-27 and nj = 18-52. Using these initial values, Table 2 is 
obtained: 








Point m’ “i n’ m’ obtained by 
I 14°8 2 1:06 9°27 Preliminary work 
2 14°5 I 0:96 18-52 Preliminary work 
3 14°6263 23°6275 1°0020 14°3239 Interpolation 
4 14°61678 23°63339 0999160 14°59087 Parallel chord 
5 14°615892 | 23°633942 0998892 14°615782 | Extrapolation 





Apparently point 5 is still on the same side of n’ = m’ as points 3 and 4 
although it is obtained from them by extrapolation. Bearing in mind the 
rapidity with which n’ varies with m’ and the closeness of the values of m’ 
and n’ for point 5, it is likely that there is asmall computational error in the 
value of n; that is unavoidable without working to more figures, so that the 
m; value does represent m, to the accuracy given. 

This example shows that when a steep branch of the curve is used the 
convergence is also very rapid. It illustrates the fact that an initial choice 
can be a good starting-point even when the revised value obtained in the 
first stage of iteration differs considerably from it. The only difficulty in 
using such branches is that they need a more careful initial search to locate 
them than one in which the revised value varies less critically with the 
starting value. 
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Without going into details, the values for m’,/, and «’ for point 5in Table 2 


lead to the following roots of the sextic equation: 
3°73570-+7 2:62628, 
4-24833 +7 0-86747, 
—2-51594+i 449266. 


These agree with the previously obtained values to within one unit in the 
last decimal place except for the second pair where the discrepancy is as 
much as four units. The real parts of the six roots in the first set add up 
correctly to —21, but it would be necessary to carry the iteration one stage 
farthey in each case to be certain of the roots to five decimal places. As in 
all iterative processes, the computations tend to become laborious as the 
need to include more and more figures increases with the degree of accuracy 


attained 


5. Conclusion 

The process described for assuring rapid convergence becomes difficult 
to apply when the value of s in (29) is close to unity and the n(m) curve lies 
nearly on top of nm = m for a considerable distance on either side of the 
value my. As may be anticipated, this difficulty arises when the sextic 
equation has several nearly equal roots or one root much larger than the 
rest that are of the same order of magnitude. 

These are the conditions in which any iterative method would be difficult 
tohandle, the reason being that the values of the roots of the sextic equation 
then depend very critically upon the precise values of the coefficients. It 
has been assumed above that even when the coefficients are given as whole 
numbers their values are exact, whereas in a practical problem their 
accuracy may depend upon the limitations of experimental measurement. 
Thus where the values of the roots are critically dependent on the values 
of the coefficients, each case must be treated on its own merits in assessing 
the accuracy to which it is justifiable to extract the roots of the equation. 

Examples of these difficult cases have been examined in which the location 
of the intersection of the n(m) curve with n = m has required a much closer 
inspection of the situation before the routine processes for convergence can 
beapplied. It is felt that the labour entailed in this more extensive explora- 
tory work does not place the method at a disadvantage compared with the 
behaviour of other iterative methods under the same difficult conditions. 

In conclusion, therefore, the method is put forward as an alternative to 
other types of iterative solution in the hope that it may prove of interest 
to those concerned with the solution of sextic equations. 





366 A NOTE ON THE SOLUTION OF THE SEXTIC EQUATION 


6. Acknowledgements 

The author is indebted to Mrs. N. A. Huttly who carried out all the 
computations, and he also wishes to acknowledge the permission of 
Marconi’s Wireless Telegraph Company Limited to publish the paper, 


REFERENCES 
. F. W. G. Otver, Phil. Trans. Roy. Soc. A, 244 (1952), 385-415. 
. Selected Government Research Reports, vol. 5 (Servo-Mechanisms). Cf. Reports 
3, 5, and 6. (H.M. Stationery Oftice, London, 1951.) 
. W. E. THomson, Wireless Engineer, 29 (1952), 256. 
. 8S. Dartineton, Bell System Technical Journal, 31 (1952), 613. 








ion of 


Paper, 











{ NOTE ON THE NUMERICAL INTEGRATION OF 
FIRST-ORDER DIFFERENTIAL EQUATIONS 


By L. FOX (The National Physical Laboratory, Teddington, Middlesex) 
Received 19 February 1953] 


SUMMARY 
This note examines methods of solving first-order differential equations with the 
ise of boundary conditions at two distinct points. The inconvenience of these 
ethods is noted, and the use is preferred of a derived second-order equation, solved 


ither by step-by-step or relaxation methods. 


1. IvreREST has been aroused by a paper of Allen and Severn (1) in which 
they remark that differential equations of the first order are not conveniently 
soluble by relaxation methods, and that it is best to use a second-order 
equation derived from the given equation of first order. The purpose of this 
ote is to add weight to that remark, and to suggest what is probably the 
best method for first-order equations, for use either by relaxation or step- 
by-step techniques. 


2. Consider the differential equation 


dy 


+S ey) = 0 (1) 


with a boundary condition specifying the value of y at some value of x. 


Suppose that this boundary point is labelled 0 in Fig. 1, and that it is 


Fic. 1. 


equired to integrate the equation by using finite-difference methods, as 
lar as the point labelled n, the constant interval between successive pivotal 
points being denoted by h. Allen and Severn quote, as the simplest approxi- 
ition to (1), the finite-difference equation 
Ys4—Yp—-1 + 2hf (2,; Y,) = 9, (2) 
but point out that one cannot solve the set of algebraic equations which 
result when equation (2) is applied at each of the points 1, 2,...,n—1. The 
ilue at the point 0 is known, but there remain n unknowns with only 
n—1) equations: if (2) is applied at the point n, another unknown, this time 
hot required, is introduced at the point n+ 1. 


[Quart. Journ. Mech. and Applied Math., Vol. VII. Pt. 3 (1954)] 
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It has been suggested to me by several workers in this field that this 
difficulty can be overcome by using formulae of type (2) at the points 
1, 2,...,2—1, and a formula at » which involves an expression for the first 
derivative using backward differences rather than central differences, and 
which therefore contains no pivotal value outside the range. The central. 
difference formula implicit in (2) is given by 

dy, i a _ 
‘ 1 1 3 1 3 1 7 } ‘ 
. dx 3(Y, +1— Yr- i)— HOY, + 39HO°Y, — 40H! O'Y, +s 5 (3 
and the corresponding backward-difference formula is 
dy, 
dx 


: ‘ 1 13 4 
h = $y,—2y,_,+4y,-2+ 4V¥y,+4V4y,4+.... (4 
Neglecting differences of orders three and higher in (4), we can replace 
(2), at the point n, by the approximate equation 
3Y —2y,, -—1 7 bY, oth f(tn, Yn) = U. (5) 

We have now n equations in n unknowns, and can obtain a solution. 

[ shall examine this method and expose its limitations by means of a 
simple example. 

3. Consider the differential equation 

dy 


y= — (6 
dx * 


with a boundary condition as yet unspecified. The difference equation 
corresponding to (2) is given by 
Yrs1—Yra— 2hy, = 9, 7) 
and this has the analytical solution 
y, = Ae™+(—1)"Be-"™, where sinh « h. (8 
To determine A and B we need two extra conditions, one of which is 
furnished by the given boundary condition. It is proposed that the second 
be given by an equation of type (5), but let us first take what would appear 
to be the best possible condition for this method and assume that we have 
been able to determine, by some unspecified method, the true value at the 
other boundary point n. We therefore take the conditions 


Yo = I, ¥, = €, (9) 
corresponding to nh = 1 and the solution e* = e”’ of the differential 
equation, and find the solution of the difference equations, in the form 


y = 4cosech(na){e"™(—e-"*-+e)+(—1)e7*(e"*—e)}  (n even)) 10) 


y = sech(na){e”*(e+e-"%)-+(—1)’e-"%(e"*—e)!  (n odd) 
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Wesee that both Ae’* and Be-'* are smooth functions, but when combined 
in the form (8) and computed at successive values of r, there will be a lack 
of smoothness caused by the factor (—1)’ multiplying Be-". 

We note further from equation (10) that the oscillation disappears com- 
pletely only when e”* = e or h = 0, and we then have the correct solution. 
It might be expected that the solution obtained by our numerical method 
is correct to the number of figures to which the function is effectively free 
from oscillation. By the latter is meant that the differences of some order 
p shall oscillate about zero with an amplitude not greater than 2?-!, the 
upper limit attainable for a polynomial ‘perturbed’ by rounding errors of 
not greater than 0-5 unit in the last figure. 

A similar analysis could be made for the case in which a condition like 
(5), instead of the second of (9), was used at the second boundary point. 
We find a similar result, this condition not affecting the form of solution 
given by {d). 

4. To illustrate, some numerical results are given below for the differ- 
ential equation (6), with finite-difference approximation (7), boundary 
conditions (9), and interval h = 0-2. Once the internal values have been 
produced, the values at external points can be calculated by step-by-step 
processes, applying (7) at the boundary and successive external points. 
These values are needed to provide central differences at internal points. 


TABLE 1 


42 
4 59° 
1*0000 632 —928 
met og 1722 
I 294 94 
2010 450 —1405 
O I°4 750 611 
155 Irss 
I SOS 547 
55 392 — 937 
957 390 
2 + 783 
2°71 )59 393 
5 I 395 
1384 


Table 1 gives the complete solution, and the oscillation in the fourth 
difference suggests that not more than two decimals are generally correct. 


lable 2 shows the functions and their differences obtained by omitting 


Bb 
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separately the terms corresponding to odd and even r, that is the functions | 
Ae’™+ Be-"® (r even) and Ae™™*— Be-"® (r odd). We see that the differences rT 
of these functions are perfectly well behaved. : 
TABLE 2 
P - y 
—o'4 +0°6674 —o6 +0°5646 
2¢ 2670 
o"o I*CO0O0 + 1600 o*2 0°8316 tr 1330 Th 
492 + 789 4CO00 + 640 pre 
tog 1°4926 2389 380 or2 1°2316 1970 
7315 1109 5970 957 waa 
oS = 22241 3558 572 06 1°8286 2927 165 do 
10873 1741 oe a 
I*2 33114 5299 ‘0 2°7183 4349 
16172 13246 se 
+ 1°6 4°9286 I°*4 +4°0429 dif 
The errors in this approximate solution, in units of the fourth decimal, 
are as follows: 39 
m 
2 0-0 0-2 0-4 0-6 0-8 1-0 
y—et 0 102 8 65 —14 0 
eq 
If, instead of the second of (9), we use, as will generally be necessary, the uy 
second boundary condition corresponding to (5), given here by to 
1-3y,—2Yy,-1+0°5y,_2 0, (11) re 
we find a solution with the following errors: ; 
x 0-0 0-2 0-4 0-6 0-8 1-0 
y—e 0 103 9 67 —13 2 
5. The possibility now arises of improving the approximation of Table | 
by applying the difference correction, that is by taking subsequent account ™ 
of the terms involving differences in equations like (3). It is at once clear ™ 
that it is not easy to decide how many differences to retain in such a formula, b 
the differences not being convergent. We are helped to some extent in that 
we require only mean odd differences, which causes some cancellation in the ) 
oscillatory terms, but the difficulty remains. For example, the difference u 
correction at the point « = 0-4 is proportional to ( 
4(301)+-3,(—247)—;1,(1032)-+..., | 
and it is not certain whether this series is really convergent. 
The difference correction itself is likely to be oscillatory, as can be seen by 
examining the mean third differences in Table 1. We can try, however, t0 
effect some improvement by incorporating it in the difference equations, 
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ind for the present illustration we will use only the third difference cor- 
rection. The errors are then reduced to the following values, expressed in 


units of the fourth decimal: 


0-0 0-2 0-4 0-6 O-S 1-0 
a 0 > 14 3 22 0 
‘he errors are smaller, but the third decimal is still inaccurate, and the 


process must be repeated several times to obtain four-figure accuracy. In 
many cases the original errors show such large oscillations that this method 
does not even converge 

If we attempt Lo improve the result corresponding to the use of (5) as the 
second boundary condition, account must be taken of third and higher 
lifterences in (4). Again we are undecided as to how many differences are 
really significant. If we use only }V*, then the ‘correction’ to y results in 
nn answer farther from the truth than the first approximation, and the 


method does not converge 


6. The correct solution will of course satisfy the complete finite-difference 


equation at an interval 0-2, and the difficulty of obtaining it is bound 
ip with the existence of a fairly large difference correction corresponding 
) the first approximation. We should get better accuracy if we could 
educe this difference correction, and for this purpose we follow Fox and 
Goodwin (2), and replace the general equation (1) by the finite-difference 


equation 
a: (q,{40” [LO ‘ oe JY, 
LhlAf(a,, y)t-f(@,—15 Up) +S (pats Ypaa)], (12) 


which is still of the three-term type, but with a much smaller difference 


rrection. For the special equation (6), the approximate form of (12) 


1— th)y,.,—(1+4A)y,_, — shy, 0. (13) 


With the boundary conditions (9) we can solve the simultaneous equations 


typified by (13), obtaining a solution whose errors, in units of the sixth 


lecima 


ire as follows. Again the oscillation is apparent, but the maximum 


error is considerably reduced. 


0-0 0-2 0-4 0-6 0-8 1-0 


If the second of (9) were not available we should have to use a formula 


similar to (5), but we should now incorporate more differences, since the 
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fifth difference is the first neglected at other points. We therefore keep 
differences up to order four in (4), obtaining in place of (11) the relation 


113 2 4 | —_ 
69 Un 4y,,- 1 9Yn-2—39n-3 T 49 n-4 0. (14) 


With (14) applied at boundary point », and (13) at all internal points, 
we find a solution with the following errors, expressed in units of the sixth 
decimal. 


x 0-0 0-2 0-4 0-6 0-8 1-0 
y—e* 0 32 14 47 36 Te 


7. We see from the foregoing that in most cases there is a restriction 
on the interval size which is a serious disadvantage of this method. It is 
well to realize that this restriction is of quite a different nature from that 
usually encountered in solving differential equations, and is caused solely 
by the presence of the factor (—1)" in equation (8). This in turn is caused 
by the use of a three-term recurrence relation for the solution of a first-orde) 
differential equation. It is, on the other hand, quite natural to use a three- 
term recurrence relation for second-order equations: no factors like (—1) 
can be present, approximate solutions do not oscillate but have convergent 
differences, and the work is no harder, the same number of equations being 
involved, each of the same three-term character. There is therefore much 
to be said in favour of the suggestion of Allen and Severn referred to in 
section 1. 

In one sense, moreover, the work will be easier. Computers know that 
relaxation methods are easier to use on second-order equations, for which 
the diagonal terms of the resultant simultaneous equations are relatively 
large, than on first-order equations, for which off-diagonal terms are usually 
larger than those on the diagonal. I am here using the term ‘relaxation’ 
in its true sense, to mean merely a certain iterative technique of solving 
simultaneous algebraic equations. 

(The different behaviour of the finite-difference equations corresponding 
to first- and second-order equations can be discussed easily for linear 
equations with constant coefficients, for which an analytical solution is 
available for the finite-difference equations. The finite-difference equation 
corresponding to the differential equation 


y' +ay = 0 
is given by Y,-44+ 2ahy,—y,_, = 9, 
and has the analytical solution 


y = Apj+ Bps, 
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where Py and P2 are roots of 
p*?+2ahp—1 = 0. 


Rither p, or p, is negative, and this causes oscillation. 
The corresponding finite-difference equation for the second-order equa- 


tion ” 
: y +ay = 9 


is given by Yr+y (2—ah?)y, Ty Q, 
and the quantities p, and p, are now roots of 
p*—(2—ah?)p+1 = 0. 


Both roots are either positive or negative, and negative only if ah? > 2. 
Thus, if a is negative, there will be no oscillation: if a is positive, oscillation 
is possible only ifhva > 2, and at this interval the finite-difference equation 
is meaningless. ) 

[ shall now add further support to the suggestion made by Allen and 
Severn by showing that one way of making the first-order equations more 
convenient for relaxation produces effectively the second-order equations. 

8. Consider the general set of equations (7). To make them convenient 
for the relaxation method we ‘normalize’ them, the equation for y, being 
obtained by multiplying each equation by the coefficient of y, contained 
therein, and adding the results. This gives the equation 


Yo + Yp-2—(2+4h?)y, = 0, (15) 


the terms in y,.,,, ¥,-; cancelling. We recognize (15) as being the approxi- 
mate finite-difference equation at an interval 2h corresponding to the 
differential equation obtained from (6) by differentiation and subsequent 
elimination of the first derivative, and given by 
d*y . 
Sy =9%. (16) 
daz 
We see now that the ‘odd’ and ‘even’ parts of the solution of the first-order 
equation have almost independent existence, satisfying separately equa- 
tions of type (16) at most points, and being connected only through the 
boundary conditions. Let us examine this connexion. 
Referring to Fig. 1, we find that if y, has the specified boundary value, 
the normalized equation for y, is given by 
y¥,—(1+4h?)y,—2hy, = 9, (17) 
which can be written in the form 
Y3g+Y_y—(2+4h?)y,+(y,—y_,—2hyy) 0, (18) 


where y_, is the value of y at the external point —1. This equation, too, we 
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recognize as being composed of the form (15) of the second-order equation 
(16), together with the form (7) of the first-order equation (6). 

If we take boundary conditions (9), so that y,, is also specified, the normal- 
ized equation for y,,_, has a form similar to (18). We conclude that, if the 
number of points in the range, including boundary points, is odd, then the 
‘even’ solution, corresponding to the points 0, 2, 4,..., satisfies the second- 
order differential equation y” = y, with specified boundary values, while 
the ‘odd’ solution, corresponding to the points 1, 3, 5,..., satisfies the same 
differential equation with conditions at the ends given by the finite- 
difference form of the first-order equation y’ = y. If the number of points 
is even, the even solution satisfies y” = y, with y specified at the beginning 
and y’ = y at the end, the odd solution satisfies y” = y with y’ = y at the 
beginning, y specified at the end. Both these solutions are perfectly smooth 
and well-behaved but they fit together smoothly only to the extent to which the 
finite-difference forms of the boundary conditions are compatible. 

If instead of the second of (9) we use a condition like (11) or (14), a similar 
analysis can be applied, but the odd and even solutions are tied together 
rather more closely, and we should expect the oscillation to be less pro- 
nounced. 

9. The above analysis gives independent arithmetical backing to the 
remark of Allen and Severn. The odd and even solutions satisfy the same 
differential equation, and will fit smoothly if we apply one condition from 
the even solution, one from the odd solution. These conditions can either 
be applied both at the same point, or one at each end of the range of integra- 
tion. For the differential equation (6), for example, with the condition given 
by the first of (8), we solve the derived equation (16) with the first of (9) used 
at the origin and the condition (6) either at the same point or at the other 
end of the range. The first method would result in an initial-value problem, 
to be solved by step-by-step methods, the second in a boundary-value 
problem, to be solved by relaxation or any technique for solving simul- 
taneous equations. 

10. I shall now indicate what seems to me the best method for applying 
these ideas to solve first-order equations. It is well known (see, for example, 
Fox and Goodwin (2)), that a second-order equation is easier to solve if the 
first derivative is absent, for we can then produce a three-term recurrence 
relation with a very small difference correction. A second-order equation 
produced from a first-order equation can always be written in this form. 
From the general equation (1), for example, we can obtain the equation 
of of 


Cy Ox 


y"—f 


0. (19) 
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snd the corresponding finite-difference equation can be written in the form 


| h2 cf | { h? of | a 5h? Of ' | 
r+ ia zy) be \"r1 et . | A9rr 6 i) og 


r J a 
h of 
10 (2) - (2) 4 (7) \_9 (20) 
») . . , 
12) Ox), OX) +1 \CX) , 1) 
where a suffix s indicates that the function is to be evaluated at x,, y,. The 


error in this equation starts with ,1,8%y,, and is therefore likely to be 


small 

If we treat the problem as one of initial value type, we need two values 
of y to start the recurrence. If the initial point has the suffix 0, we have y, 
but need to be able to determine y, or y_,. This we can do by applying the 
boundary condition’ (1) at the point 0. Because the differential equation 
has a small difference correction, we must use an accurate expression for (1). 
This we achieve by using at the point 0 the equation (12), which has a small 
orrection starting with ,446°y), and involves only three adjacent pivotal 
values. By applying (20) also at 0, we have two equations for the determina- 
tion of y_, and y,, and the step-by-step method can proceed. 

If we treat the problem as one of boundary-value type, we use (20) at 
ill internal points, and (20) and (12) applied at the second boundary point 


to provide a workable boundary condition. 


11. Applied to linear differential equations, this method is quick and 
accurate, and, if necessary, difference corrections can be applied with con- 
fidence. For equation (6), for example, with boundary condition given by 
the first of (9), we can solve by boundary-value techniques from the following 


simultaneous equations 





6:10y(0-2)+-2-99y(0-4) = —2-99 
29y(C 2)—6-10y(0-4) + 2-99y(0-6) 0 
2-99y (0-4) —6-10y(0-6)+ 2-99y(0-8) 0 
2-99y(0-6)—6-10y(0-8)+-2-99y(1-0) 0 
2-99y(0-8)—6-10y(1-0)+ 2-99y(1-2) 0 
3°20y(0-8)—0-80y(1-0)+ 2-80y(1-2) 0 

(21) 


[he error in this solution, in units of the sixth decimal, is as follows: 


0-0 0-2 0-4 0-6 0-8 1-0 
18 29 42 


e 


57 
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This is an improvement on the corresponding result of section 6 as regards 
the magnitude of the error, but of greater importance is the fact that the 
latter is now a very smooth function, devoid of oscillation, and the difference 
correction can clearly be applied with confidence. I give below the table 
of differences of this approximation, and the calculated values of the 
difference corrections which should be added to the left-hand sides of equa- 
tions (21). The quantities labelled A (= 51,5°y) should be added to the first 
five equations of (21) and the quantity A, (= 3 d5y) to the last of (21), 


x y 10°A 
—c'4 0°67¢302 
+148420 
—o2 0818722 4+- 32858 
181278 + 7276 
0'0 ~=—- T'00C000 40134 +1609 
221412 8885 +361 
+02 = 1°221412 49019 1970 + Fo +03 
270431 10855 431 
o4 =: 1"491843 59874 2401 Ic3  +0°4 
330305 13256 534 
c6 =: 1822148 73130 2935 116 +05 
403435 16191 650 
o8 2°225583 89321 3585 144 06 
492756 19776 794 
10) =. 2718339 109097 4379 + 174 +07 (108A, = 9°8) 
601853 24155 +968 
I°2 3320192 133252 + 5347 
735105 + 29502 
1°4 4°C55297 + 162754 
+897559 


16 47953156 
The resulting corrections to y in units of the sixth decimal are as follows: 


a 0-0 0-2 0-4 0-6 0-8 1:0 
dy 0 —10 — 20 —32 — 46 —62 


The difference correction corresponding to dy is clearly negligible, and we 
should have the correct answer. The remaining error in the sixth decimal 
is caused by the ill-conditioning of the linear equations: the correction to 
y is at least six times the major difference correction A,, and to obtain the 
sixth decimal more accurately we should need to work throughout with at 
least one extra figure. 


12. We saw in the last paragraph that the first-order equation can be 
solved conveniently by deriving a second-order equation, using the given 
condition at one boundary and a suitable approximation to the first-order 
equation at the same point or at the other end of the range. The 
first method I described as an initial-value technique, the second as 4 
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houndary-value technique: they might be distinguished by the terms 
‘ecurrence’ and ‘relaxation’. It may legitimately be queried whether 
relaxation has any point in its favour, since recurrence is invariably easier 
to perform. In most problems relaxation has no advantage, but in certain 
cases it will give better precision in the computed pivotal values. 
a] 
Consider, for example, a linear differential equation of the form 
y' +yf(x) = g(x). (22) 
The solution consists of a particular integral P(x) and a complementary 
function C(x), and . P 
y P(x)+A C(2), (23) 
where A is an arbitrary constant whose value is determined from the 
boundary condition on y. Suppose that C(x) is a rapidly increasing function 
fx, that P(a) varies much more slowly, and that the conditions of the 
problem are such that in the required solution A = 0, so that no multiple 
of the complementary function is present. When we perform the recurrence 
method, we shall round off our pivotal values to a certain number of figures, 
with the result that a small multiple of C(x) is introduced and the error of 
our solution increases with increasing x. By relaxation, on the other hand, 
we are solving simultaneous equations and applying a boundary condition 
it each end, and the error increase may be avoided. 
As an example, consider the solution of the differential equation and 
boundary condition given by 
y' —12y+1le* 0), y(0) is (24) 
Deriving the second-order equation, we find that the general equation 
corresponding to (20) has here the form 
Ipay(1— 12h?)+-y,_,(1 — 12h?) —y,(2+ 120h?) 
143},2 Ir Ir r— oF 
—M3h?(10er+-err+1+e%—1) (25) 
ind that corresponding to (12) is given by 
Yray(1 4h) y,-4(1-44 4h)— l6hy, -- Ih (4e? + efr+it gtr—-1), (26) 
Applying the recurrence method at an interval h = 0-2, we find the 
following values of y: 
0-0 0-2 0-4 0-6 0-8 1-0 
j 11-0000 1:2214 11-4918 1-8218 2-2214 2-6644 
The relaxation method gives the following result: 
1-0000 ‘2214 1:-4918 1-8220 2-2255 2-7184 
and for comparison the corresponding four-decimal values of e* are 
appended: 


10000 11-2214 14918 1-8221 2-2255  2-7183 
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I am indebted to Professor D. R. Hartree and Dr. E. T. Goodwin for 
useful discussions on parts of this paper, which has been written as part 


of the research programme of the National Physical Laboratory, and is 


published with the permission of the Director. 
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TABLES OF THE INTEGRALS AND DIFFERENTIAL 
COEFFICIENTS OF Gi(+a2) AND Hi(—z2) 


By M. ROTHMAN (The Northern Polytechnic, Holloway, N.7) 
Received 7 May 1953] 


1. Introduction 


Ir was shown in a previous paper (1) that the solution of the problem of 
in infinite plate under a uniform inclined loading depended on tables of the 
integrals of the functions Ai(+ 2), Bi(-+2), and Gi(+2) and the former 
two integrals were tabulated for x 0(0-1)10. 

Scorer (2) has tabulated Gi(+2) and the associated function Hi(—z), 
which can be interpolated more easily for negative values of the argument. 
In this paper the integrals and the derivatives of Gi(+ a) and Hi(—2) are 
given to seven decimal places for x 0(0-1)10. 

It can also be shown that 


1 T 


Gi(x) = 4Bi(x)+Ai(x) | Bi(u) du — Bi(x) [ Ai(u) du 


0 0 
: “4 (1.1) 
Gi’ (2) LBi'(2)+-Ai'(x) | Bi(u) du Bi'(x) | Ai(u) du 
0 0 
Gi’ (x)—ax Gi(z) 1/7 


where accents denote derivatives, so that tables of Gi(x) and Gi'(a2) form 
in alternative to those of | Ai(x) dx and { Bi(x) dz. The functions 
Gi'(+a2) and Hi'(—x) have therefore been included in this paper, 
the functions Gi(+2) and Hi(—2) being given by Scorer (2). Though the 
tables given only extend to x 10, for uniformity with (1) and (2) the 
author has computed values to x 20, which are available if required. 


2. Construction of the given tables 
The integrals were obtained by numerical integration of 10 place tables 


of Gi(+2z) and Hi(—z2), using the formula 
\ | ] (x) da w| (fo + ti) ae of, t Jot 54f,—...]. (2.1) 


rhe first differences were added up to give successive values of the 
integrals with increasing 2. The results were then checked by using the 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 3 (1954)] 
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integrated form of the asymptotic expansions of Gi(+ x) and Hi(—z) 
which give the required degree of accuracy from 2 = 8-0 onwards, i.e, 


Deo, 2). 8 8! 
Gi(7) ~ +( Benen Seas Specie Sa 
wie 2 ae" See 


a7 


Hi(—2x) ~ -(: —— ats a | 


wie 2 Sa" 3.62%" 


i9¢ 
(2.2) 


and a final check was obtained by comparing several results in the table 
with those computed from the formulae: 


4 2y+In3 1 2! 5! \ | 
| de (inx—2 | ad i] 
p vo. c” 
0 


[ Hi(—é) dé = a > 4. “(ine 5 | i -) 
é = 
where y is Euler’s constant. 
The tables of Gi'(+ax) and Hi’(—x) were obtained by numerically 
differentiating the 10 place tables of Gi(+) and Hi(—~2) using the formula 
why = hafy— tH O%fo + sor Ofo—---s (2.4) 
and the results obtained were compared with those obtained by differen- 
tiating the asymptotic expansions (2.2) where possible. 
I wish to express my thanks to J. C. P. Miller for his helpful advice and 
to F. W. J. Olver of the National Physical Laboratory for the formulae (2.3). 
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M. ROTHMAN 
dGi(+ x) 1Gi(+ x) 

dx 8, x dx 82, 

0-0 0°14942 94 205 95 50 0°01 355 95 4 35 
I 11866 92 231 11 I 01325 19 } 02 
2 09020 4I 247 59 2 01268 45 3 69 
3 06420 05 2560 43 3 O1215 40 3 41 
4 04074 94 255 57 4 o1165 76 3 09 
05 0°01987 36 25512 | §°5 O’OI11G 21 2 87 
6 00154 04 247 06 6 ©1075 53 2 64 
7 01432 59 235 34 7 01034 49 2 43 
8 02784 97 220 g2 8 00995 88 223 
9 03910 50 204 4° 9 00959 50 207 
1-0 0°04843 78 18686 60 0°00925 19 1 89 
I 05554 32 105 57 I 00592 77 175 
2 OO150 31 150 20 2 00502 13 I 62 
3 00575 05 132 10 3 00533 11 I 51 
4 06867 52 11477 4 o0805 60 I 41 
15 0°07042 06 829 65 0°00779 50 I 31 
6 o7118 11 82 g2 6 00754 71 I 22 
Zz Vill o!1 65 52 7 00731 14 Bf 
8 07034 86 55 97 8 00708 69 107 
9 o0fg02 4 44 50 9 00057 31 99 
2-0 fe 34 31 70 0°00666 g2 g2 
I 25 44 I 00047 45 oo 
2 17 76 2 00625 86 SI 
3 11 17 3 Ooo! od 77 
4 571 4 00594 07 72 
2°5 0°05452 02 L312/| 75 0°00577 75 635 
6 05208 22 2 56 6 00562 17 64 
7 04930 54 5 49 7 00547 20 60 
8 04670 83 7 80 8 00532 83 5% 
9 O4412 51 9 49 9 ©0519 04 54 
30 0°04163 59 107 8-0 000505 79 51 
I 03925 31 In 55 I 00493 05 45 
2 03698 51 1200 2 00450 79 17 
3 03453 06 1219 3 00409 0O 43 
4 03250 90 12 16 4 00457 04 {I 
35 0°0 3090 41 11 94 8-5 0°00446 69 jo 
6 o2g11 78 11 63 6 00436 14 35 
7 02744 7° 111s 7 00425 97 35 
& 2588S S« 10 6S ~ 00416 15 25 
9 02443 65 10 05 9 00406 65 31 
4°0 fe 51 g0 00397 52 32 
I gl I 00355 65 30 
2 22 2 00350 I4 27 
3 774 3 00371 57 25 
4 716 4 00363 55 25 
45 Oo 664 95 0°00356 14 25 
6 611 6 00345 65 23 
7 01596 10 5 04 7 00341 39 24 
8 + 01521 00 5 19 8 00334 37 21 
9 OI451 og 477 9 00327 50 20 

0°00 320 
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LES OF 


0-0 . yO Oe 50 0°7 52904 35 12 06 
I 205 50 I 75916 06 11 65 
2 242 34 2 760526 09 II 21 
3 219 15 3 77124 9! 10 52 
4 198 61 4 77712 91 10 45 

0°5 © 42 55 0°75290 46 10 10 
6 164 22 6 78857 91 975 
7 149 01 7 79415 OF 9 44 
8 136 97 8 79963 87 913 
9 5 47 9 80502 00 8 82 

1-0 I5 15 6:0 0°51033 30 55 
I 05 93 I Sisscsos 5 30 
2 7 59 2 52005 50 5 
3 0 12 3 57392 7 
4 4 S3071 55 7 
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Theoretical Elasticity 


By A. E. GREEN and W. ZERNA 50s. net 


The first part of this book is concerned with the general theory of finite elastic 
deformations. Special emphasis is placed on the solution of problems for incom- 
pressible isotropic bodies using a general strain-energy function. Classical 
infinitesimal theory, deduced as a special case of general theory, is used in the 
rest of the book, which mainly deals with plates and shells. Complex variable 
methods are used in the two-dimensional theory and applications for both ee 
and aeolotropic bodies. The later chapters contain a general theory of thin shells, 
a theory of shallow shells and membrane shells, and general methods of solution 
of problems for cylindrical shells and shells of revolution. Tensor and vector 
notations, used consistently throughout the book whenever appropriate, are sum- 
marized in Chapter 1. 
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By J. A. SCHOUTEN Second edition 30s. net 


‘This subject is expounded in a monumental manner. . . . The treatment is full 
and exhaustive. . . . There is probably no other book which covers the same 
ground in such extent and in such detail.’ NATURE, reviewing the first edition 
published in 1951. 
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